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ABSTRACT 

Three  problems  in  the  stability  of  viscoelastic  flows  have 
been  theoretically  investigated. 

The  case  of  plane  Couette  flow  has  been  solved  by  the  classi- 
cal energy  method  of  Orr.  For  the  method  to  be  applicable  it  has 
been  found  that  the  concept  of  elastic  potential  energy  of  a visco- 
elastic liquid  has  to  be  introduced,  and  an  expression  has  been 
found  for  it.  For  two-dimensional  disturbances  of  any  magnitude 
the  presence  of  elasticity  has  been  found  to  stabilize  the  flow. 
This  result  suggests  that  the  viscous  sublayer  thickness  must 
increase  during  the  so-called  Toms  phenomenon.  It  has  been  found 
that  the  results  are  identical  for  three  fluids,  namely  the  second- 
order  fluid  and  the  Walters  fluids  A*  and  b',  although  the  normal 


stress  behavior  of  these  fluids  are  quite  different.  This  insensi- 
tivity to  the  constitutive  equation  has  been  shown  to  be  a result 
of  the  two-dimensionality  of  the  problem. 

Linear  stability  theory  has  been  used  to  investigate  the 
stability  of  two  flows  of  the  Oldroyd  fluid,  namely  the  plane 
Poiseuille  flow  and  the  circular  Couette  flow.  For  both  flows  the 
presence  of  elasticity  has  been  found  to  destabilize  the  flow  and 
to  increase  the  critical  wavenumber. 

As  an  aside,  it  has  also  been  shown  that  the  sign,  and  magni- 
tude, of  the  Weissenberg  effect  may  depend  on  the  method  of 
observation,  namely  through  the  shape  of  the  free  surface  or 
through  the  height  difference  of  radial  tappings. 
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CHAPTER  1 
INTRODUCTION 


1.1  The  Toms  Phenomenon 

Great  interest  was  created  in  the  mechanics  of  non-Newtonian 
fluids  after  the  discovery  of  the  so  called  Toms  phenomenon.  Toms 
(1948)  discovered  that  when  small  quantities  of  long  chain  linear 
polymers  having  molecular  weights  of  the  order  of  five  or  six  figures 
are  added  to  a turbulent  shear  flow,  substantial  reduction  (by  factors 
of  two  or  three)  in  the  pressure  drop  and  drag  are  observed  at  the 
same  flow  rate.  There  was  some  conjecture  that  the  reason  might  be 
that  the  addition  of  the  polymer  actually  stabilized  the  flow,  trans- 
forming the  turbulent  flow  into  a laminar  one,  resulting  in  a reduction 
in  skin  friction  at  the  same  Reynolds  number. 

Many  theoretical  and  experimental  attempts  were  subsequently  made 
to  investigate  the  stability  characteristics  of  non-Newtonian  fluids. 

It  was  found  that  no  great  change  in  the  critical  Reynolds  number 
occurs  due  to  the  non-Newtonian  properties  of  the  fluid.  Therefore, 
the  conclusion  was  that  the  Toms  phenomenon  cannot  be  due  to  the 
delayed  transition  of  a non-Newtonian  fluid  [see,  for  example,  Fong  § 
Walters  (1965),  Tlapa  5 Bernstein  (1970)].  Moreover,  it  has  been 
observed  that  the  flow  does  actually  remain  turbulent  during  the  Toms 
phenomenon. 

No  single  explanation  of  the  Toms  phenomenon  has  been  universally 


accepted  among  the  scientific  community.  It  was  shown  by  Lumley  (1964) 
that  the  purely  viscous  non-Newtonian  fluid  (Reiner-Rivlin  fluid) 


r 
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cannot  cause  drag  reduction;  he  concluded  that  dependence  of  the 
stress  on  the  history  of  the  strain  rate  is  necessary,  in  other  words 
a viscoelastic  material  is  necessary.  For  speculations  as  to  the 
mechanism  responsible  for  the  Toms  phenomenon, see,  for  example, 

Lumley  (1967,  1969),  and  Fabula,  Lumley  § Taylor  (1966). 

Although  the  conjecture  that  the  Toms  phenomenon  is  due  to 
delayed  transition  is  now  known  to  be  false,  the  investigation  of  the 
instability  of  a viscoelastic  fluid  by  the  so  called  "energy  method" 
can  nevertheless  be  connected  to  the  phenomenon.  For  example  the 
experimentally  observed  thickening  of  the  viscous  sublayer  during  the 
Toms  phenomenon  will  be  explained  in  the  present  work  by  means  of  a 
stability  analysis. 

1.2  Notes  on  the  Two  Methods  of  Stability  Analysis 

The  subject  of  stability  deals  with  the  question  of  the  ability  or 
the  inability  of  a physical  system  to  sustain  itself  against  perturba- 
tions to  which  it  is  subject.  In  the  area  of  fluid  mechanics  it  is 
known  that  laminar  flows  occur  at  low  Reynolds  numbers  and  turbulent 
flows  at  higher  Reynolds  numbers.  Thus,  a typical  problem  of  hydro- 
dynamic  stability  asks  the  question:  At  what  Reynolds  number  does  the 

laminar  flow  become  unstable  so  that  disturbances  can  grow  (and 
eventually  break  it  up  into  disorderly  turbulent  motion)? 

There  are  two  different  methods  for  investigating  the  stability 
of  a disturbed  flow  - the  method  of  small  disturbances  (linear  theory) 
and  the  method  of  finite  disturbances  (nonlinear  theory).  In  the 
method  of  small  disturbances  one  linearizes  the  resulting  equations  by 
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neglecting  terms  quadratic  or  higher  in  the  disturbances  and  their 
derivatives.  This  method  has  so  far  received  more  attention.  For 


example  among  the  major  treatises  on  hydrodynamic  stability, 
Chandrasekhar  (1961)  deals  exclusively  with  this  approach,  and  Lin 
(1955)  deals  very  largely  with  this  approach.  Although  quite  success- 
ful in  a large  number  of  cases  "the  linear  theory  suffers  from  the 
defect  that  one  cannot,  in  principle,  make  judgements  regarding  the 
growth  potential  of  finite  disturbances.  Thus,  one  cannot  say  for 
certain  that  a given  flow  will  remain  stable  if  disturbed  under  con- 
ditions judged  favorable  by  the  linear  theory"  [Joseph  (1966)]. 

To  see  this  point  more  clearly,  let  us  suppose  that  it  is  possible 
to  decompose  the  motion  into  elementary  components  so  that  the  flow  can 
be  considered  as  a set  of  interacting  elementary  nonlinear  oscillators 
[see  Monin  § Yaglom  (1971),  page  8S] . "In  each  of  these  oscillators, 
due  to  the  influx  of  energy,  self-excited  oscillations  may  arise.  The 
possibility  of  such  oscillations  arising  is  determined  by  the  relation- 
ship between  the  energy  E+  acquired  by  the  oscillator  and  the  energy  E~ 
lost  by  it  for  various  amplitudes  a of  the  oscillations  (Fig.  1.1). 

If  E > E+  for  all  amplitudes  (Fig.  1.1a)  then  the  oscillations  will 
evidently  be  damped  for  any  initial  amplitude,  and  the  system  will  be 
stable  for  any  disturbances.  If  E-  < E+  for  a1  < a < aQ,  but  E~  > E+ 
for  a < ax  or  a > aQ  (Fig.  1.1b),  then  the  oscillations  with  initial 
amplitude  a < a^  will  be  damped,  but  those  with  initial  amplitude 
a > will  increase,  until  their  amplitude  attains  the  equilibrium 
value  aQ.  In  this  case,  the  system  will  be  stable  to  small  distur- 
bances but  unstable  to  disturbances  of  sufficiently  large  amplitude  ... 


Figure  1.1  Different  variants  of  the  dependence  of  the  energy  gained 
and  lost  by  an  oscillator  on  the  amplitude  of  the  oscilla- 
tion [after  Monin  § Yaglom  (1971)]. 
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Finally,  if  E*  > E”  for  any  amplitude  below  a , however  snail,  the 

0 

system  will  be  unstable  to  infinitely  small  disturbances  (i.e., 
absolutely  unstable)  and  will  practically  always  be  in  the  regime  of 
self-excited  oscillation  with  amplitude  a 

0 

There  are  other  deficiencies  in  the  linear  theory.  It  is  quite 
clear  that  the  loss  of  stability  does  not  in  itself  constitute  a 
transition  to  turbulence,  and  the  linear  theory  can  at  best  describe 
only  the  very  beginning  of  the  process  of  transition.  Moreover,  for 
certain  important  flows,  for  example  the  Couette  flow  between  two 


parallel  planes  or  Poiseuille  flow  in  a tube,  the  small  disturbance 
theory  predicts  that  the  flow  is  stable  to  all  infinitesimal  distur- 
bances at  all  Reynolds  numbers.  Also,  in  the  case  of  a Poiseuille  flow 
between  two  parallel  planes  the  linear  theory  gives  a critical  Reynolds 
number  which  is  considerably  higher  than  that  at  which  transition 
actually  occurs  in  a real  flow. 

A nonlinear  theory  of  the  initiation  of  turbulence,  which  takes 
into  account  the  finite  size  of  the  disturbances,  has  received  a lot  of 
attention  during  the  last  fifteen  years.  For  a recent  review  article 
on  this  topic  see  Stuart  (1971);  see  also  the  books  of  Betchov  § 
Criminale  (1967),  Monin  5 Yaglom  (1971)  and  the  earlier  review  article 
on  general  hydrodynamic  stability  by  Stuart  (1963). 

In  addition  to  the  nonlinear  theory  of  the  type  just  mentioned, 
there  is  another  simple  approach  to  the  investigation  of  the  stability 
of  a flow  with  respect  to  disturbances  of  finite  amplitude.  This  is 
the  so  called  energy  method,  which  is  really  classical  and  originated  in 
the  writings  of  Osborne  Reynolds  (1895).  He  supposed  that  a disturbance 
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motion  of  finite  amplitude  is  already  in  existence,  and  derived  the 
following  disturbance  energy  equation  for  a Newtonian  fluid  (for 
derivation,  see  Section  3.1): 


_d_ 

dt 


2 


/ ^ dV  - - / 

V V 


u.u.U.  . dV  - v 
i J i,3 


/ 

V 


u . . u . . dV 
i,J  i,3 


(1.1) 


where  u^  is  the  disturbance  velocity,  IL  is  the  basic  flow  velocity,  v 
is  the  kinematic  viscosity,  V is  the  volume  of  integration,  and  a comma 
denotes  differentiation  with  respect  to  space  coordinate.  The  left 
hand  side  of  equation  (1.1)  represents  the  rate  of  increase  of  distur- 
bance kinetic  energy  throughout  the  volume  V;  the  first  term  on  the 
right  hand  side  represents  the  rate  of  extraction  (or  production)  of 
disturbance  energy  through  the  interaction  of  the  Reynolds  stresses 

u.u.  with  the  gradient  U.  . of  the  basic  flow;  the  last  term  in  (1.1) 

* 3 1 » J 

is  always  negative  and  represents  viscous  dissipation. 

Reynolds  criterion  is  that  the  flow  is  stable  or  unstable  depending 
on  whether  the  kinetic  energy  of  the  disturbance  motion  will  decrease  or 
increase  with  time,  i.e., 

* 

- / u.u.U.  .dV  < v / u.  .u.  .dV  (for  stability) 

V 1 3 y 1 1 J 1 > J 

(Production)  (Dissipation) 

The  equality  sign  in  the  above  relationship  corresponds  to  "neutral 
stability,"  for  which  the  rate  of  change  of  disturbance  kinetic  energy 
is  zero. 

The  first  careful  use  of  the  energy  method  was  made  by  Orr  (1907). 


He  considered  all  disturbances  which  give  rise  to  instantaneous 
stationarity  of  disturbance  kinetic  energy,  for  which  (1.1)  reduces  to 


- I Vjui,jdv  * “ I ui.jui,/v  »•« 

For  a given  function  u^x),  a value  of  v can  be  found  from  (1.2).  The 

question  raised  by  Orr  is  the  following:  "What  function  u^  makes  v 

the  maximum?"  For  a larger  value  of  viscosity  the  motion  must  be 

stable.  For  a Newtonian  plane  Couette  flow  having  two-dimensional 

2 

disturbances  he  found  a critical  Reynolds  number  of  T6  /v  = 44.3,  where 
T is  the  shear  rate,  25  is  the  distance  between  the  plates,  and  v is 
the  kinematic  viscosity.  For  plane  Poiseuille  flow  with  centerline 
velocity  he  obtained  a critical  Reynolds  number  Umc/v  - 87. 

In  a modem  formulation  of  the  energy  method,  Serrin  (1959)  sought 
to  find  conditions  under  which  the  kinetic  energy  of  the  disturbance 
motion  will  surely  decrease;  that  is,  he  found  sufficient  conditions 
for  stability  under  arbitrary  three-dimensional  disturbances.  Serrin's 
method  depends  on  finding  the  upper  or  lower  bounds  on  the  various 
term.-i  in  the  disturbance  energy  equation,  and  the  limits  of  stability 
given  by  estimates  of  this  kind  do  not  depend  on  the  specific  details 
of  the  geometrical  configuration  or  on  the  distribution  of  the  basic 
field  variables.  This  generality  is  obtained  at  the  expense  of  rather 
more  conservative  estimates  than  can  be  obtained  by  using  all  of  the 
available  information. 

Joseph  (1966)  applied  the  variational  method  of  Orr  to  a heated 

2 

plane  Couette  flow  and  obtained  the  stability  boundary  of  R + Ra  = 

2 

1708,  where  R = 2T5  /v  and  Ra  is  the  Rayleigh  number.  For  zero 

2 

Rayleigh  number  limit,  the  criterion  reduces  to  T5  /v  = / 1 708/ 2 = 20.6. 
This  value  is  about  one  half  of  that  given  by  Orr;  the  reason  for  the 
discrepancy  is  that,  whereas  Orr  assumed  that  the  two-dimensional 
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disturbances  are  more  dangerous  than  three-dimensional  ones,  Joseph 
considered  longitudinal  roller-type  disturbances  and  proved  that  they 
are  more  unstable. 

The  energy  method  has  been  subject  to  extensive  criticism  for 
quite  some  time.  For  example,  Lin  (1955),  Serrin  (1959),  Joseph  (1966) 
and  Schlichting  (1968)  state  that,  in  this  method,  one  establishes 
stability  relative  to  arbitrary  disturbances  while,  in  reality,  only 
those  satisfying  the  hydrodynamical  equations  need  be  considered. 

Lumley  (1971b)  in  a recent  comment  on  the  energy  method  has  shown  that 
this  criticism  is  incorrect,  "since  functions  are  admitted  for  review 
at  a fixed  time  only;  any  velocity  field  at  a fixed  time  satisfies  the 
hydrodynamical  equations,  which  serve  to  determine  its  evolution,  i.e., 
the  time  derivative."  In  the  same  paper  Lumley  also  shows  that  the 
disturbance  motion  in  the  energy  method  must  necessarily  be  nonlinear, 
and  that  the  disturbance  structure  must  change  because  the  local  growth 
rate  must  be  different  from  point  to  point. 

The  one  particular  criticism  about  the  energy  method,  that  it  gives 
low  values  of  the  critical  Reynolds  number,  is  true.  But  Lumley  (1971b) 
comments  that  this  is  to  be  expected  since  the  class  of  disturbances 
considered  is  larger  than  that  in  the  linear  theory.  Another  (perhaps 
the  primary)  reason  for  the  lower  value  is  the  fact  that  in  this  method 
one  does  not  exclude  "the  possibility  that  disturbances  which  have,  at 
the  instant  of  analysis,  positive  growth  rates,  may  ultimately  decay 
(Fig.  1.2)  ...  Thus,  one  says  that  at  a given  Reynolds  number  the 
energy  method  can  ascertain  stability,  but  cannot  determine  instability." 


GLOBAL 
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Figure  1.2  Initial  growth  of  a disturbance  which  ultimately  decays 
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A great  success  of  the  energy  method  is  that  it  gives  a nonzero 
value  of  the  critical  Reynolds  number  for  the  plane  Couette  flow 
problem,  whereas  the  small  disturbance  theory  predicts  this  flow  to  be 
stable  at  all  Reynolds  numbers  [see,  for  example,  Petrov  (1938)]. 

Recent  work  of  Joseph  5 Hung  (1971),  however,  suggests  that  even  the 
roller-type  disturbances  in  parallel  flow,  shown  by  the  energy  method 
to  be  initially  unstable,  are  ultimately  stable. 

1. 3 Survey  of  Theoretical  Work  Done  on  Viscoelastic  Fluid  Stability 

Considerable  work  has  been  done  on  the  stability  of  viscoelastic 
fluids  during  the  last  fifteen  years,  and  a general  survey  is  given 
below.  (For  an  explanation  of  terms  peculiar  to  viscoelastic  fluids, 
like  "second-order  fluid,"  "Walters  fluid,"  or  "Oldroyd  fluid"  etc., 
reference  may  be  made  to  Chapter  2.) 

(i)  Plane  Couette  Flow:  For  a plane  Couette  flow  of  a Newtonian 

fluid  the  linearized  stability  equations  admit  only  the  trivial  (that 
is,  zero)  solution  for  the  disturbance  quantities  at  any  Reynolds 
numbers,  implying  that  the  flow  is  always  stable.  Giesekus  (1966) 
observed  that  this  is  not  true  in  presence  of  viscoelastic  parameters, 
and  he  obtained  a stability  criterion  for  the  plane  Couette  flow  of  a 
viscoelastic  fluid  using  the  linear  stability  theory.  We  however  feel 
that  the  application  of  the  linear  theory  to  the  plane  Couette  flow  may 
be  misleading,  even  of  a viscoelastic  flow,  and  a thorough  discussion 
of  this  point  is  given  in  Section  1.5. 
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(ii)  Horizontal  Layer  of  Fluid  Heated  from  Below  (Benard 

Stability) : When  a static  layer  of  Newtonian  liquid  is  heated  from 

below,  instability  sets  in  at  a critical  value  of  the  so-called 
Rayleigh  number.  It  is  found  that  the  principle  of  exchange  of 
stabilities  holds  for  this  case,  and  the  neutral  state  is  characterized 
by  steady  cellular  convection.  Intuitively,  one  may  guess  that  this 
may  not  be  true  in  a viscoelastic  liquid  in  which  the  elastic-like 
behavior  of  the  fluid  may  give  rise  to  "overshooting"  and  an  oscilla- 
tory type  of  behavior  at  the  marginal  state  ("overstability").  Green 
(1968)  considered  this  possibility  for  a linearized  Oldroyd  fluid,  and 
Vest  5 Arpaci  (1969)  examined  it  for  a linearized  Maxwell  fluid.  They 
both  came  to  the  conclusion  that  overstability  is  possible  in  a visco- 
elastic fluid,  but  that  for  ordinary  thicknesses  of  common  liquids  the 
relaxation  times  would  have  to  be  tremendously  high  for  overstability 
to  take  place;  therefore,  a simple  experiment  to  demonstrate  the  over- 
stability is  not  feasible  with  currently  available  fluids.  They  also 
found  that  the  elasticity  of  the  fluid  has  a destabilizing  effect  on  a 
liquid  heated  from  below.  This  is  true  both  in  the  sense  that  oscilla- 
tory convection  can  occur  at  a lower  critical  Rayleigh  number  than  does 
stationary  convection,  and  that  this  critical  number  decreases  as  the 
viscoelastic  parameter  increases. 

(iii)  Plane  Couette  Flow  With  Heating:  Herbert  (1963)  considered 

a very  general  Oldroyd  fluid  (which  displays  variation  of  viscosity 
with  shear  rate)  and  found  that  the  Rayleigh  number  is  decreased  for 
all  nonzero  rates  of  strain;  this  destabilization  was  attributed  solely 
due  to  the  variation  of  the  apparent  viscosity  with  shear  rate. 
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Mclntire  5 Schowalter's  (1970)  analysis  is  valid  for  both  a 

"generalized  second-order  fluid"  and  an  integral  equation  which  has  its 

basis  in  a molecular  model  and  accounts  for  network  junctions  rupturing 

at  a certain  critical  rate;  they  found  that  the  flow  is  stabilized  or 

destabilized  according  as  the  second  normal  stress  difference  a = 

2 

S22  - S33  is  negative  or  positive. 

(iv)  Circular  Couette  Flow:  This  is  the  stability  of  flow 

between  two  concentric  rotating  cylinders.  Jain  (1957)  considered  a 
Reiner-Rivlin  (or  Stokesian)  fluid,  but  errors  in  Jain's  work  were 
pointed  out  by  Graebel  (1961)  who  found  that,  for  fluids  with  a positive 
coefficient  of  cross-viscosity,  the  critical  Taylor  number  can  be 
appreciably  smaller  than  for  the  corresponding  flow  of  a Newtonian 
fluid.  Thomas  fj  Walters  (1964a)  found  a decrease  in  the  critical 
Taylor  number  due  to  the  non-Newtonian  parameter  for  the  Walters  liquid 
b'  with  short  memory,  whereas  Fong  (1965)  found  that  the  opposite  is 
true  for  the  Walters  liquid  a'  (without  short  memory) . Thomas  5 Walters 
(1964b)  also  found  a destabilization  for  a Maxwell  fluid  B.  Datta 
(1964,  1966),  Rao. (1964)  and  Ginn  5 Denn  (1969)  considered  a second- 
order  fluid;  the  general  conclusion  was  that  the  second  normal  stress 
difference  a (proportional  to  B + 2y)  and  the  diameter  to  gap  ratio 
are  crucial:  a2  > 0 always  results  in  destabilization,  but  c>2  < 0 
gives  rise  to  a destabilization  in  wide  gap  and  stabilization  in  narrow 
gap. 

The  possibility  of  overstability  of  a circular  Couette  flow  was 
considered  by  Beard,  Davies  5 Walters  (1966) . They  considered  a 
Maxwell  fluid  B of  stress  relaxation  time  \ , viscosity  u,  density  p 
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and  gap-width  6.  They  concluded  that  for  values  of  K = yA^/pi  up  to 
K = 0.13,  the  only  Tc  (critical  Taylor  number)  which  exists  corresponds 
to  a stationary  stability  mode.  In  this  region,  the  presence  of 
elasticity  in  the  liquid  is  a significant  destabilizing  agent.  For 
higher  values  of  K,  a Tc  exists  for  both  a stationary  and  an  over- 
stable mode.  Since  the  critical  Taylor  number  associated  with  the 
overstable  mode  was  always  lower  than  that  associated  with  the 
stationary  mode,  they  concluded  that  overstability  sets  in  at  K = 0.13 
and  is  present  for  all  higher  values  of  K.  The  critical  Taylor  number 
associated  with  the  overstable  mode  was  found  to  decrease  steadily 
with  increasing  K,  so  that  "highly  elastic  Maxwell  liquids  can  be  very 
unstable  indeed." 

(v)  Flow  in  a Curved  Channel:  This  is  the  flow  in  a curved 

channel  due  to  a pressure  gradient  acting  along  the  channel.  Thomas 
Walters  (1965)  studied  the  stability  of  this  flow  for  a Walters  B 
fluid  (without  assuming  short  memory),  but  the  results  were  given  for 

a Maxwell  fluid  B which  is  a special  case  of  the  fluid  they  considered. 
They  found  a destabilization,  together  with  an  increase  in  the  critical 
wavenumber  of  the  convective  cells,  due  to  the  elasticity  of  the  liquid. 

(vi)  Plane  Poiseuille  Flow:  This  is  the  stability  of  flow  driven 

by  a pressure  gradient  between  two  parallel  planes.  Fong  § Walters 
(1965)  treated  it  for  a Walters  b'  fluid  with  short  memory,  Chun  5 
Schwarz  (1968)  for  a second-order  fluid,  Tlapa  § Bernstein  (1968,  1970) 
for  a generalized  Maxwell  fluid  (errors  in  Tlapa  5 Bernstein's  work 
will  be  pointed  out  in  Section  4.11),  and  Jones  5 Walters  (1968)  for  a 


third-order  fluid  which  exhibits  a variation  of  viscosity  with  shear 
rate.  All  of  them  concluded  a destabilization  due  to  the  non-Newtonian 
parameters . 


(vii)  Dangers  of  Using  Second-Order  Fluid  Model:  Quite  a few 

articles  have  pointed  out  the  dangers  of  using  the  second-order  con- 
stitutive equation  in  a linear  stability  analysis.  Coleman,  Duffin  5 
Mizel  (1965)  have  shown  that  all  solutions  for  unsteady  parallel  flows 
of  a second-order  fluid  are  unstable.  Pipkin  (1966)  commented  that 
this  apparent  instability  is  "an  interesting  absurdity  which  arises 
when  the  ...  [slow  flow]  approximation  is  treated  as  if  it  were  exact." 
Gupta  (1967)  showed  that  a layer  of  second-order  fluid  at  rest  between 
plane  boundaries  is  in  unstable  equilibrium;  Craik  (1968)  also 
reiterated  the  same  fact.  Platten  5 Schechter  (1970)  showed  that  a 
uniform  rigid  translation  (U  = constant  everywhere)  of  second-order 
fluid  is  unstable  at  all  Reynolds  numbers.  These  results  contradict 
the  earlier  mentioned  works  of  Fong  5 Walters  (1965)  and  Chun  5 Schwarz 
(1968)  who  calculated  stability  envelopes  and  did  not  mention  any  such 
trouble  with  second-order  fluids.  The  paradox  was  explained  by 
Mclntire  (1971).  He  considered  a plane  Couette  flow  with  superposed 
temperature  gradient,  and  found  curves  of  marginal  stability  as  shown 
in  Fig.  1.3. 

Note  carefully  that  for  high  wavenumbers  the  flow  is  unstable  at 
all  Rayleigh  numbers.  If  Fong  5 Walters  (1965)  and  Chun  5 Schwarz 
(1968)  had  examined  the  very  high  wavenumber  perturbations,  they  would 
have  found  an  instability  at  all  Reynolds  numbers  also.  However,  it 
would  be  quite  incorrect  to  conclude  that  a real  viscoelastic  liquid  is 


Figure  1.3  Stability  diagram  of  heated  Couette  flow  of  second-order 
fluid,  for  given  values  of  Reynolds  number  and  non- 
Newtonian  parameters  [after  Mclntire  (1971)]. 
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always  unstable  at  high  wavenumbers,  because  the  second-order  fluid  is 
not  applicable  near  this  high  wavenumber  stability  boundary.  The 
reason  is  that,  as  noted  in  Fig.  1.3,  the  growth  rate  is  continuous 
through  zero  across  the  ordinary  (low  wavenumber)  stability  envelope, 
whereas  for  the  high  wavenumber  instability  the  growth  rate  is  dis- 
continuous and  very  large.  In  order  for  the  second-order  model  to  be 
applicable  the  growth  rate  of  the  disturbance  must  be  small  compared  to 
a "natural  time"  of  the  fluid. 

I 

(viii)  Application  of  Nonlinear  Methods:  In  all  the  works  cited 

' 

so  far  the  linear  theory  of  stability  analysis  was  applied.  Feinberg  5 
Schowalter  (1969a,  1969b)  seem  to  be  the  first  to  consider  disturbances 
of  finite  magnitude  and  apply  the  energy  method  to  the  investigation  of 
stability  of  viscoelastic  flows.  They  followed  the  technique  of  Serrin 
and  determined  sufficient  conditions  for  stability  of  fluid  motions 
constitutionally  described  by  the  infinitesimal  theory  of  visco- 
elasticity. Like  Serrin 's  work  on  Newtonian  fluids,  their  results 
were  also  overly  conservative. 

Denn,  Sun  § Rushton  (1971)  carried  out  a finite  amplitude  analysis 
of  the  Secondary  motion  in  a circular  Couette  flow,  using  the  energy 
equilibrium  approach  of  Landau  and  Stuart.  Assuming  a second-order 
fluid  model,  they  calculated  that  the  extra  torque  following  the 
instability  is  less  than  that  for  the  Newtonian  fluid  at  the  same 
relative  (supercritical)  Taylor  number.  They  suggest  that  this  might 
somehow  explain  the  drag  reduction  in  turbulent  pipe  flow  of  polymer 


solutions . 
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1.4  Survey  of  Experimental  V.'ork  Done  on  Viscoelastic  Stability 


All  the  experimental  work  on  the  stability  of  polymer  solutions 

have  been  performed  on  the  circular  Couette  flow  set-up,  with  the  outer 

cylinder  stationary.  The  results  have  been  quoted  in  terms  of  the  non- 

2 3 2 

dimensional  Taylor  number  T = 6 R /v  where  U and  R are 

11  1 ) 

respectively  the  angular  velocity  and  the  radius  of  the  inner  cylinder, 
<5  is  the  gap  width,  and  v is  the  kinematic  viscosity.  The  results  of 
the  various  observers  are,  however,  sharply  contradictory. 

For  example  Giesekus  (1966)  found  an  extremely  strong  destabiliza- 
tion. For  a 4%  solution  of  polyisobutylene  in  decalin  he  found 
instabilities  in  which  the  rotational  speed  was  about  600  times  lower 
(Taylor  number  approximately  3 x 105  times  smaller.')  than  in  a New- 
tonian liquid.  For  a 1%  solution  of  Al-naphthenate  in  decalin  he  found 
the  critical  rotational  speed  to  be  about  30  times  lower  (Taylor  n amber 
about  900  times  smaller)  compared  to  a Newtonian  liquid.  Giesekus  also 
reported  that  the  critical  wavenamber  of  the  convective  cells  increases 
(that  is,  the  spacing  decreases)  due  to  the  addition  of  polymers. 

Miller  (1967)  also  found  a definite  destabilization,  but  reported  that 
the  critical  wavenumber  decreases  with  the  addition  of  polymers. 

On  the  other  hand  Karlsson  5 Griem  (1966)  and  Rubin  5 Elata  (1966) 
found  definite  stabilization  even  for  dilute  polymer  solutions  for 
which  the  deviation  from  Newtonian  behavior  is  otherwise  small.  Rubin 
G Elata  also  reported  that  the  cellular  spacing  remains  relatively  con- 
stant at  its  Newtonian  value. 

Merrill,  Mickley  5 Ram  (1962)  found  that,  if  v is  taken  as  the 
viscosity  at  the  critical  shear  rate  then  the  solution  of  polymer 
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macromolecules  does  not  show  any  significant  change  in  the  critical 
Taylor  number  compared  to  micromolecular  Newtonian  liquids.  However, 
they  also  commented  that  if  v were  taken  as  the  viscosity  at  zero  shear 
rate,  then  their  results  would  indicate  a substantial  destabilization. 

Denn  5 Roisman  (1969)  found  that  for  dilute  polymer  solutions  the 
critical  Taylor  number  may  be  more  or  less,  depending  on  the  diameter 
ratio  of  the  cylinders. 

1.5  Critical  Examination  of  Giesikus'  Work  on  Plane  Couette  Flow 
Instabi li  1L 

It  was  mentioned  in  Section  1.3  that  Giesekus  (1966)  analyzed  the 
stability  of  a viscoelastic  plane  Couette  flow  by  means  of  linear 
stability  theory.  A critical  examination  of  this  work  is  given  below. 

Giesekus  considered  plane  Couette  flow  as  a limiting  case  of 
circular  Couette  flow  when  the  radii  of  the  walls  go  to  infinity.  That 
is,  he  considered  that  the  disturbances  are  periodic  in  the  x^-direction 
and  are  in  the  form  of  "longitudinal  rollers"  invariant  in  the  x^- 
direction  (x^  is  taken  in  the  direction  of  the  flow  and  x2  is  taken 
perpendicular  to  the  flow).  The  stability  equation  obtained  was 

222  2 22  2 * * 
t(D  -k  ) + k y(D  -k  ) - k Y ]u»  = 0 (1.3) 

* A 

with  the  boundary  conditions  i|i  = Di>  = 0 at  x^  = 0 and  1.  Here  all  the 
variables  are  dimensionless,  D = d/dx2>  k is  the  wavenumber  of  the  dis- 
turbance, $ is  the  amplitude  of  the  stream  function  in  the  x^x^-planc 
and  y = - G(F-G)r2/2p2  and  y*  = - Gp62r  /p2  are  two  elastic  parameters. 
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T is  the  shear  rate  of  the  basic  flow,  u is  the  viscosity,  p is  the 

density,  6 is  the  gap  width,  and  F and  G are  two  normal  stress 

functions  defined  to  be  related  to  the  normal  stresses  in  the  fluid  by 

2 2 

the  equations  = - p + T F,  = - p-r  G,  and  S = - p,  where  p is 
the  pressure. 

For  a Newtonian  fluid  y = y =0,  and  the  only  solution  of  (1.3) 
which  satisfies  the  boundary  conditions  is  the  trivial  solution  $ = 0. 
This  means  that  for  a Newtonian  fluid  the  transverse  yz-flow  is  always 
stable  at  any  Reynolds  numbers.  Giesekus  observed  that  this  is  not 
true  with  y and  y*  present,  and  he  solved  (1.3)  and  obtained  the 
following  stability  criterion:  for  a highly  elastic  liquid  the 

instability  is  to  be  expected  if  y - 1,  whereas  for  slightly  elastic 
liquids  the  instability  is  to  be  expected  if  y*  = 68.  (It  of  course 
follows  that  G will  have  to  be  negative  to  cause  instability.) 

It  is  apparent  that  Giesekus  confined  his  attention  only  to  the 
transverse  motion.  However,  the  longitudinal  disturbance  velocity  may 
be  crucial  in  the  study  of  stability.  For  example,  for  longitudinal 
roller  type  disturbances  in  a parallel  flow  of  a Newtonian  fluid, 

Joseph  § Hung  (1971)  have  shown  that  the  transverse  motion  monotonically 
decays,  while  the  longitudinal  motion  (Uj)  may  initially  increase  if  the 
Reynolds  number  is  high  enough.  Lumley,  in  a private  communication,  has 
shown  that  the  same  fact  is  also  true  of  linear  disturbances  in  a 
Newtonian  fluid.  This  is  explained  below. 

The  linearized  equations  of  motion  for  the  disturbance  in  this 


case  are 


I 

't 


3u  2 

——  + Tu  = V 7 u 
dt  2 1 
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(1.4) 


(1.5) 


2 2 

where  7=3  /3x^3x^,  and  is  the  stream  function  for  the  transverse 
motion.  The  boundary  conditions  are  ^ = \p  = dtp/dx^  at  = i 5 , 
Multiplying  Eq.  (1.5)  by  ip,  integrating  over  the  flow  volume,  and  using 
the  divergence  theorem  and  boundary  conditions  many  times,  results 


1 

2 


_d_ 

dt 


/ 


2 

q dV  = - v 


/ 


u>  dV 


(1.6) 


2 2 2 2 
where  q = u + u is  the  transverse  energy  and  w = 7 ij>  is  the 

vorticity.  This  shows  that  the  transverse  energy  must  monotonically 

decay,  which  of  course  must  happen  since  there  is  no  source  term  in 

Eq.  (1.5). 

To  investigate  the  ^-motion,  one  must  realize  that  any  initial 
distributions  of  u1  and  ip,  consistent  with  boundary  conditions,  are 
possible  - the  equations  will  determine  their  evolution.  For  example, 
suppose  a ^ is  found  that  satisfies  the  boundary  conditions;  then 
initially  let  us  choose  u « - a|i  where  c is  a constant,  which  also 
satisfies  the  boundary  conditions  on  u^  After  multiplying  Eq.  (1.4) 
by  Uj,  integrating  over  the  flow  volume,  and  using  the  divergence 


theorem  and  boundary  conditions  we  get 
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_d_ 

dt 


2 
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CN 

f 2 

* 

3ul 

u dV  = - 

u dV  - v 

1 c 

1 

3x. 

l 

dV 


(1.7) 


It  is  apparent  that  although  the  second  term  on  the  right  side  is 
always  negative,  for  large  enough  r the  growth  rate  may  be  made 
positive  initially. 

The  longitudinal  motion  having  thus  been  shown  to  be  crucial,  the 
work  of  Giesekus  becomes  questionable.  It  seems  that  a method  which 
deals  with  all  three  components  of  velocity  in  a single  equation,  and 
treats  simultaneously  all  parts  of  the  flow,  would  be  preferable  in 
this  case.  The  classical  energy  method  is  such  a one;  it  has  the 
additional  advantage  that  it  is  unnecessary  to  assume  the  disturbances 
small  at  the  beginning. 


1.6  Estimating  Viscous  Sublayer  Thickness  Through  Energy  Method 

G.  I.  Taylor  (1916)  suggested  using  the  energy-method  type  of 
analysis  for  the  purpose  of  estimating  the  thickness  of  the  viscous 
sublayer  next  to  the  wall  in  a turbulent  flow,  and  Lumley  (1971b) 
recently  applied  the  idea  to  estimate  the  thickness  of  the  viscous  sub- 
layer in  a Newtonian  fluid.  Their  argument  is  as  follows:  a viscous 

sublayer*  exists  along  the  walls  of  a turbulent  flow  because  the 
proximity  of  the  wall  makes  the  flow  here  dominated  by  viscosity.  Now 
imagine  a turbulent  plane  Couette  flow  (Fig.  1.4);  the  two  viscous 


Note  that,  strictly  speaking,  it  is  incorrect  to  speak  of  a "laminar" 
sublayer  since  it  consists  of  unsteady  disturbances,  although 
dominated  by  viscosity. 
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sublayers  in  this  flow  are  expected  to  have  the  same  structure  as  in 
any  other  turbulent  flow.  If  the  viscosity  of  the  fluid  is  gradually 
increased,  other  things  being  equal,  the  sublayer  disturbances  will  not 
be  able  to  sustain  themselves  without  increasing  in  thickness.  The 
sublayers  will  therefore  go  on  increasing  in  thickness  until  they  are 
face  to  face.  This  is  the  so-called  "critical  state"  of  stability 
studies,  at  which  point  a further  increase  in  viscosity  will  damp  out 
all  disturbances  and  produce  a steady  laminar  flow. 

These  considerations  can  be  phrased  as  a well  defined  extremization 
problem:  consider  the  class  of  large  disturbances  having  instantan- 

eously stationary  total  energy;  what  is  the  equation  obeyed  by  that 
member  which  can  exist  with  the  largest  viscosity  (other  things  being 
equal) ? 

The  procedure  of  estimating  the  viscous  sublayer  thickness  is  now- 
clear:  at  the  neutral  or  critical  state  the  viscous  sublayers  are  face 

to  face,  and  half  the  distance  between  the  plates  is  the  sublayer 
thickness.  Using  Orr's  solution  (which  considered  only  two-dimensional 
disturbances)  for  the  critical  Reynolds  number  and  noting  that  by 
definition  u^  = /vF  , one  gets  <5+  = Su^/v  = /FT^/v  = /44.5  = 6.65; 
using  Joseph's  solution  (which  considered  longitudinal  roller-type 
disturbances)  one  gets  6+  = /20 . 6 = 4.54.  These  values  are  in  fairly 
good  agreement  with  experimental  values  of  a Newtonian  viscous  sublayer. 

In  the  present  work  these  same  ideas  will  be  used  in  an  effort  to 
explain  the  thickening  of  the  viscous  sublayer  due  to  viscoelasticity, 
as  observed  in  the  Toms  phenomenon. 
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From  the  discussion  of  Section  1.5  it  will  be  seen  that  the  linear 
stability  analysis  is  not  applicable  to  plane  Couette  flow,  and  a major 
objective  of  the  present  study  is  to  obtain  a stability  criterion  for 
plane  Couette  flow  by  means  of  the  energy  method.  Since  the  critical 
Reynolds  number  of  this  flow  gives  a method  of  estimating  the  viscous 
sublayer  thickness  of  any  turbulent  flow  (as  discussed  in  Section  1.6), 
one  objective  of  this  study  would  be  to  predict  whether  the  presence  of 
viscoelasticity  should  increase  or  decrease  the  viscous  sublayer  thick- 
ness. Of  course  it  is  known  from  the  Toms  phenomenon  that  the  sublayer 
thickness  increases  after  the  addition  of  polymers  to  the  solution. 

From  the  survey  of  Section  1.3  it  will  be  seen  that  the  Oldroyd 
constitutive  equation,  which  is  probably  the  most  reliable  of  all  the 
simple  constitutive  equations,  has  not  been  applied  to  the  investiga- 
tion of  the  linear  instability  of  plane  Poiseuille  flow  and  circular 
Couette  flow.  Results  of  course  exist  for  the  Maxwell  fluid  for  both 
of  these  flows,  but  the  Maxwell  fluid  does  not  exhibit  the  experi- 
mentally observed  effect  of  strain  rate  relaxation.  Besides,  we 
detected  some  errors  in  Tlapa  5 Bernstein's  (196S,  1970)  solution  of 
the  plane  Poiseuille  problem,  which  will  be  pointed  out  in  Section  4.11. 
One  objective  of  the  present  study  will  therefore  be  to  apply  the 
Oldroyd  constitutive  equation  to  the  investigation  of  the  instability 
of  the  plane  Poiseuille  flow  and  the  circular  Couette  flow. 

The  energy  considerations  for  the  various  viscoelastic  fluid 
models  do  not  seem  to  have  received  any  attention  in  the  literature. 
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One  objective  of  the  present  study  will  therefore  be  to  find 
expressions  for  the  rate  of  change  of  elastic  potential  energy  and 
the  rate  of  viscous  dissipation  for  the  commonly  used  fluid  models 
such,  as  the  second-order  fluid,  Maxwell  fluid,  Walters  fluid  and 
Oldroyd  fluid.  It  will  be  seen  in  Chapter  3 that  these  expressions 
will  be  necessary  in  the  application  of  the  energy  method  to  plane 
Couette  flow  of  a viscoelastic  fluid. 

As  an  aside,  some  unsteady  flow  problems  will  be  solved  for  the 
Oldroyd  fluid;  these  do  not  seem  to  have  been  solved  in  the  literature 
so  far. 


CHAPTER  2 


VISCOELASTIC  CONSTITUTIVE  EQUATIONS  AND  THEIR  ENERGY  CONSIDERATIONS 


2.1  Non-Newtonian  Liquids 


An  incompressible  Newtonian  fluid  is  one  for  which  the  stress 


tensor  T„  is  given  by 


T.  . = - PS. . + S. . 

iJ  ij  ij 


(2. 


with 


S..  = 2uE.. 


(2.2) 


where  P is  an  indeterminate  hydrostatic  pressure,  5^  is  the  Kronecker 


delta,  S. . is  called  the  extra-stress  tensor,  E. . = -i (U . . ♦ U.  .)  is 

ij  i,J  3,iJ 


the  rate  of  strain  tensor,  and  p is  a constant  (viscosity)  which  can 
only  be  a function  of  temperature.  Any  liquid  whose  stress-strain  rate 
relation  is  different  from  the  linear  relation  (2.2)  is  called  a non- 
Newtonian  liquid. 

For  some  non-Newtonian  liquids  it  may  be  possible  to  explain  its 
properties  completely  in  terms  of  a single  viscosity  coefficient,  if 
this  is  considered  to  be  a function  of  the  local  rate  of  strain  of  the 
material : 

S = n(r)r  (2.3)  ' 


where  S is  the  shear  stress,  r is  the  local  rate  of  shear  and  n is  an 
even  function,  called  the  (variable)  coefficient  of  viscosity  which  may 
be  a function  of  the  rate  of  shear.  A liquid  whose  behavior  at  any 
point  where  the  motion  is  one  of  simple  shearing  (not  necessarily 
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steady)  conforms  with  (2.3)  may  be  called  a non-Newtonian  viscous 
liquid.  For  these  liquids  the  mode  of  flow  under  different  conditions 
can  be  worked  out  without  great  difficulty  in  principle.  But  there  are 
several  rheological  phenomena  observable  in  liquids  which  cannot  be 
explained  in  terms  of  this  simple  theory;  in  particular  the  effect  of 
elasticity  of  shape  of  the  liquid  becomes  important.  These  are  called 
viscoelastic  liquids,  and  polymer  solutions,  molten  plastics,  biological 
fluids,  protein  solutions  (e.g.,  egg  white),  and  most  industrial  liquids 
in  the  form  of  suspensions  and  emulsions  are  viscoelastic.  Most  visco- 
elastic liquids  also  display  a variable  viscosity  with  shear  rate. 


2.2  Some  Properties  of  Viscoelastic  Liquids 

Some  important  properties  of  viscoelastic  liquids  are  described 
below.  For  a more  complete  description  with  photographs,  see  the  film 
notes  of  Markovitz  (1968). 

(a)  Variable  Viscosity:  For  most  viscoelastic  liquids  the  shear 

stress-shear  rate  relation  is  nonlinear,  implying  a coefficient  of 
viscosity  which  varies  with  the  shear  rate.  With  some  liquids  such  as 
some  polymer  solutions  the  viscosity  decreases  with  increasing  shear 
rate;  these  are  called  shear-thinning  or  pseudoplastic.  The  decreasing 
viscosity  is  thought  of  as  arising  from  a change  in  structure  of  the 
material.  With  some  suspensions,  the  viscosity  increases  with  shear 
rate;  this  is  called  shear- thickening  or  dilatant  behavior. 

(b)  Normal  Stress  Effects:  In  steady  simple  shearing  flow  between 


infinite  parallel  plates,  the  normal  stresses  on  an  infinitesimal 
volume  element  in  the  three  orthogonal  directions  are  not  equal  for  a 
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viscoelastic  liquid,  i.e.,  t T22  f ^ ; but  in  a Newtonian  liquid, 

they  are  equal.  The  same  is  also  true  for  the  normal  stresses  in  the 
flow  between  rotating  cylinders  or  in  the  torsional  shear  flow  between 
rotating  parallel  disks.  As  a result,  the  distribution  of  normal 
stresses  can  be  quite  different  for  a viscoelastic  liquid  from  what  it 
is  for  the  Newtonian  fluid. 

For  example,  in  circular  Couette  flow  in  the  annulus  between  an 
outer  cylinder  and  a coaxial  rotating  shaft,  a viscoelastic  liquid 
tends  to  climb  up  the  inner  cylinder;  this  is  called  the  (positive) 
Weissenberg  effect  [heissenberg  (1947)],  which  is  in  contrast  with  the 
Newtonian  case  in  which  the  free  surface  is  parabolic.  It  can  be 
shown  [see  Truesdell  § Noll  (1965),  p.  455]  that  the  positive  IVeissen- 
berg  effect  is  possible  only  if  the  normal  stresses  at  a point  are  not 
all  equal. 

(c)  Expansion  of  Jet:  A jet  of  viscoelastic  liquid  swells  as  it 

emerges  from  a conduit;  this  is  called  Merrington's  effect.  It  can  be 
shown  that  [see,  for  example,  Jaunzemis  (1967),  p.  538]  for  this  to 
happen  T22  > T , where  1 is  the  direction  of  flow. 

(d)  Stress  Relaxation:  In  a perfectly  viscous  liquid  the  stresses 

will  become  hydrostatic  as  soon  as  the  flow  is  stopped.  In  a visco- 
elastic liquid  the  stress  returns  to  its  hydrostatic  value  gradually . 

(e)  Elastic  Recovery:  If  the  stress  is  reduced  to  its  hydrostatic 

value,  the  fluid  flow  will  immediately  stop  in  a perfectly  viscous 
liquid;  but  in  a viscoelastic  fluid  it  will  do  so  gradually . 

(f)  Other  Elastic  or  Memory  Effects:  In  a purely  viscous  liquid 

the  current  value  of  the  stress  is  determined  from  a knowledge  of  the 
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current  state  of  deformation  alone.  In  a viscoelastic  liquid,  the 
stress  is  a (nonlinear)  function  of  the  deformation  which  the  material 
has  experienced  previously,  i.e.,  is  a function  of  the  deformation 
gradient  history' . As  a consequence  of  this  memory  effect,  such  liquids 
show  elastic-like  behavior.  For  example,  a "silicone-putty  ball 
resting  on  the  ground  flows  into  a puddle  even  under  such  a small  force 
as  that  of  gravity,  if  given  sufficient  time.  But  during  the  short 
time  of  a rapid  impact,  it  behaves  like  an  elastic  solid  because  it 
does  not  have  a chance  to  "forget"  its  previous  spherical  shape.  The 
significant  relaxation  times  for  stresses  in  silicone  putty  are  much 
less  than  the  time  of  rapid  impact.  Thus  the  putty  may  exhibit  a range 
of  behavior  from  elastic  to  viscous,  depending  on  how  the  character- 
istic time  of  the  experiment  compares  with  the  significant  relaxation 
times"  [Markovitz  (1968)]. 

2.3  Simple  Fluid 

In  Sections  2.3  through  2.7,  some  important  constitutive  equations 
will  be  discussed.  We  shall  always  limit  ourselves  to  incompressible 
fluids,  whether  or  not  the  incompressibility  condition  is  mentioned 
explicitly.  We  shall  start  with  the  very  general  concept  of  simple 
fluid. 

The  concept  of  simple  fluid  was  originated  by  Noll  (1958),  and 
developed  subsequently  by  Coleman  and  Noll  [see,  for  example,  Coleman  8 
Noll  (1961)].  They  defined  a simple  material  as  a substance  for  which 
the  stress  is  determined  by  a knowledge  of  the  entire  history  of  the 
deformation  gradient  (strain).  A simple  material  is  called  a simple 


fluid  if  it  has  the  additional  property  that  all  local  states  are 
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equivalent  in  response,  with  all  observable  differences  in  response 
being  due  to  definite  differences  in  history;  that  is,  it  has  the 
property  that  every  configuration  is  an  undistorted  configuration. 

Consider  a material  point  P which  occupies  the  position  x at  the 
present  time  t.  Let  us  suppose  that  at  an  earlier  time  t,  the  particle 
occupied  the  position  For  the  dependence  of  £ on  x,  t and  r we 
write 

S = 

The  function  X may  be  called  the  displacement  function  relative  to  the 
configuration  at  time  t.  The  spatial  gradient  of  Xt(x,x)  with  respect 
to  x determines  the  relative  deformation  gradient  tensor  at  time  x 
relative  to  the  configuration  at  time  t: 

S,M  * 5x?t<x'T’ 

We  now  define  a symmetric  tensor  through  the  relation 

ct(x)  = - l 

where  the  superscript  T denotes  "transpose"  and  I denotes  the  unit 
tensor.  G^fx)  is  a measure  of  the  "strain"  of  the  configuration  at 
time  x relative  to  the  configuration  at  time  t;  for  all  x therefore 
) may  be  looked  upon  as  the  entire  history  of  the  deformation  at 
the  material  particle  P. 

It  is  convenient  to  regard  Gt(x)  as  a function  of  the  time  lapse 
s between  x and  t,  i.e.,  s = t-x  where  o < s < °°.  Noll  showed  that  for 
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simple  fluids  the  extra  stress  S is  related  to  the  history  of  the 
motion  through  an  equation  of  the  form 

00 

S(t)  = F [Gt(s)]  (2.4) 

s=o 

where  F is  a tensor-valued  functional,  i.e.,  an  operator  (which  need 
not  be  linear)  mapping  an  entire  tensor-valued  function  G into 

at 

another  tensor  S. 

Simple  fluids  can  exhibit  such  phenomena  as  stress-relaxation, 
non-Newtonian  viscosity,  and  normal  stress  effects.  It  may  be  remarked 
that  a better  name  for  simple  fluids  would  probably  be  general  fluid, 
because  they  include  perfect  fluids,  Newtonian  fluids,  the  Reiner- 
Rivlin  fluids,  the  Rivlin-Erickscn  fluids,  and  the  fluids  based  or.  the 
theory  of  finite  linear  viscoelasticity  as  special  cases. 

In  recent  years  much  has  been  learned  about  the  mechanical 
behavior  of  general  incompressible  simple  fluids.  It  has  been  found 
that  there  are  several  steady  flow  problems  which  can  be  solved  for 
simple  fluids  in  general,  using  no  special  constitutive  assumptions 
beyond  the  definition  (2.4);  simple  shearing  flow,  Poiseuille  flow  and 
circular  Couette  flow  are  examples. 

2.4  Rivlin-Ericksen  Fluid 

Rivlin  6 Ericksen  (1955)  considered  the  theory  of  isotropic 
materials  for  which  the  extra-stress  depends  on  only  the  spatial 
gradients  of  velocity,  acceleration,  second  acceleration,  ...  , (N-l)th 
acceleration.  They  showed  that  for  an  incompressible  material  of  this 


r 
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type,  the  extra  stress  must  be  given  by  an  isotropic  function  of 


the  tensors  M = 1,  2,  ...  , N,  i.e.,  an  N-th  order  Rivlin- 


Ericksen  fluid  is  given  by 


S . . *-  h . . (A  , A , ...  , A*,. ..n) 

ij  x3  (l)fcd  (2)k£*  (N)k^ 


(2.5) 


where  the  A,.,,...  called  the  Rivlin-Ericksm  tensors,  are  given  by 
(N  j ij  — 


A a U.  . + U.  . 

Cl)  ij  *.J  J*1 


(2.6) 


A - %N-l)ij  + n + A U 

(N)ij  Dt  (N-l)ik K,j  (N-l)jk 


Here  D/Dt  = 3/3t  + U^C  ) It  can  be  shown  [see  for  example  the 
recent  review  article  of  Rivlin  § Sawyers  (1971)]  that  the  Rivlin- 


Ericksen  constitutive  equation  is  an  approximation  to  Noll's  constitu- 


tive equation  for  an  incompressible  isotropic  simple  fluid  with  memory. 


provided  that  the  Cauchy  strain  Gt(t)  can  he  expanded  as  a Taylor 


series  about  the  time  t at  which  the  stress  is  measured.  Many  important 


constitutive  equations  are  special  cases  of  (2.5).  A perfect  fluid 


corresponds  to  the  case  = 0,  and  for  a Newtonian  is  assumed  to 


depend  linearly  on  A^^  only,  i.e.,  •-*  where  u is  a 


constant.  In  case  of  the  Reiner-Rivlin  fluid,  the  extra  stress  depends 


on  only  A^j^,  but  is  taken  as  a general  (nonlinear)  function  of  it; 
for  that  fluid  equation  (2.5)  and  the  representation  theorem  for 


isotropic  tensor  functions  of  one  symmetric  rensor  variable  yield 


J 
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S..  = f A,  ...  + f A,  . . , A , . , . 

1J  1 (OlJ  2 (Uik  (l)kj 


(2.7) 


where  f and  f are  scalar  functions  of  the  invariants  of  A,,....  The 

1 2 (l)lj 

Reiner-Rivlin  constitutive  equation  (2.7)  is  purely  viscous  and  has  no 
elastic  properties. 

Although  the  general  Rivlin-Ericksen  fluid  (2.5)  accounts  for 

shear- dependent  viscosity  and  the  normal  stress  effects,  it  has  the 

serious  shortcoming  of  not  accounting  for  gradual  stress  relaxation. 

When  A......  = 0 for  M = 1,  2.  ...  , N,  the  extra  stress  S..  of  a Rivlin- 

(M)ij  ij 

Ericksen  fluid  cannot  change  in  time,  contrary  to  the  behavior  of  real 
viscoelastic  fluids  in  stress-relaxation  experiments.  Stated 
differently,  "when  a Rivlin-Ericksen  fluid  is  brought  to  rest  suddenly 
its  memory  does  not  fade  gradually,  but  rather  precipitously." 

For  steady  viscometric  flows  the  Rivlin-Ericksen  fluid  simplifies 
considerably.  For  such  flows  only  the  first  two  tensors  A.  ...  and 

(l)ij 

A,,....  are  different  from  zero.  Thus  if  no  other  tests  beyond  those 
(2)lJ  7 

resting  upon  steady  viscometric  flows  were  available,  we  should  be 

unable  to  distinguish  between  the  general  Rivlin-Ericksen  fluid  and  the 

special  one  in  which  S. . = h. . (A,  ..  n.  A,  ») . If  we  assume  that  h. . 

ij  ij  (l)  k-t  (2)k  l1  ij 

is  a polynomial  function,  by  means  of  the  representation  theorem  of 
isotropic  matrix  polynomial  of  two  matrices  it  can  be  shown  that 


2 2 2 2 
S=aA  + a A + a A +aA  +a(AA  +AA)+a(A,A  + A A ) 
= 1 = 1 2 = 2 3=1  W = 2 5-1  = 2 =2=1  6 =1=2  =2=l' 


2 2 2 2 2 2 
+ a (A  A + A A ) + a (A  A + A A ) 

7 =1=2  =2=1  8 =1=2  =2=1 


■ 


34 

where  a , ...»  a are  polynomials  in  the  ten  invariants  formed  from 

A and  A„ . 

.1  =2 

2.5  Second-Order  Fluid 

Due  to  the  generality  of  the  simple  fluid,  it  is  not  surprising 
that  not  many  flow  problems  can  be  solved  with  it.  Approximations  to 
the  simple  fluid  then  becomes  necessary. 

To  the  concept  of  simple  fluid  Coleman  5 N'oll  (1960)  added  the 
notion  of  gradually  fading  memory,  i.e.,  the  idea  that  deformations 
which  occurred  in  the  distant  past  should  have  less  effect  on  the 
present  value  of  the  stress  than  the  deformations  which  occurred  in  the 
recent  past.  For  this  they  constructed,  for  an  arbitrary  history 
G(  ),  a new  history  J(  ) for  each  a , 0 < <x  < 1 , as  follows 

G ^ (s)  = G(cts)  0 < s < °° 


(Note  that  we  have  dropped  the  subscript  t on  G , which  should  cause  no 

r s 

confusion.)  Gl  } ( ) is  therefore  the  same  history  as  G(  ),  but  carried 
out  at  a slower  rate;  it  can  therefore  be  called  the  retardation  of 
G(  ) by  the  factor  a.  Coleman  5 Noll  (1960)  have  shown  that  the 
stresses  corresponding  to  G^(  ) obey  the  approximate  formulae 

S(a)  = ♦ 0(cO  (2.8) 

« = U J 


;(“)  - I,  a (®)  + ftA(a)  ,(a)  .(a) 

: ' “hu  %)  ho  vh2) 


+ 0(0 


(2.9) 
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where  a|^  = a^A^>  A^  are  the  Rivlin-Ericksen  tensors,  and  p,8  and 
Y are  material  constants.  Equation  (2.8)  tells  us  that,  for  a simple 
fluid  with  fading  memory,  the  Newtonian  fluid  furnishes  a complete 
first-order  correction,  for  viscoelastic  effects,  to  the  theory  of 
perfect  fluids  in  the  limit  of  slow  motion.  Equation  (2.9)  exhibits  all 
second-order  corrections  for  viscoelastic  effects.  The  equation 
obtained  by  striking  out  the  error  term  0(a3)  in  (2.9)  is 


] 


UA 


(l)ij 


eA0)ikA(l)kj 


yA 


(2)ij 


(2.10) 


This  is  called  a second-order  fluid,  and  can  be  looked  upon  as  an 
approximation  to  the  simple  fluid. in  slow  motions. 

Note  that  the  theory  of  Reiner-Ri vlin  fluids  (2.7)  would  not  yield 
the  term  in  (2.9)  involving  y;  hence  that  theory  cannot  give  a complete 
second-order  correction  to  the  Newtonian  fluid. 

The  main  practical  interest  in  second-order  fluids  lies  in  its 
simplicity,  and  quite  a number  of  hydrodynamic  problems  can  be  solved 
with  it.  In  several  situations  where  the  theory  of  Newtonian  fluids 
gives  rise  to  linear  partial  differential  equations  for  the  velocity 
field  so  also  does  the  theory  of  second-order  fluids.  Reference  may  be 
made  to  Markovitz  8 Coleman  (1964)  for  a good  summary  of  the  derivation 
and  applications  of  the  second-order  constitutive  equation.  It  is 
shown  there  that  the  fluid  displays  the  normal  stress  effect  and  the 
positive  Weissenberg  effect.  However,  the  second-order  fluid  does  not 
display  the  gradual  stress  relaxation  and  the  shear-dependent  viscosity. 


Truesdell  (1965)  showed  that  the  second-order  fluid  is 
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indistinguishable  bv  visconetric  experiments  from  a fluid  with  convected 
elasticity.  In  the  theory  of  fluids  with  convected  elasticity  the 
stress  is  regarded  as  a function  of  the  strain  of  the  present  con- 
figuration with  respect  to  one  occupied  by  the  fluid  at  a certain  fixed 

* * 

time  t before  the  present  time.  The  time  t , which  is  a material 

property,  may  be  called  the  response  time  of  such  a fluid.  Truesdell 
showed  that  in  viscometric  flows  a fluid  of  second  order  is  indistin- 
guishable from  a particular  fluid  of  convected  elasticity  with  response 

time 

' t*  = - (2.11) 

u 

Since  u > 0,  this  time  is  positive  if  and  only  if  y < 0.  This  fact 
affords  us  a provisional  interpretation  for  the  coefficient  y.  The 
second-order  fluid  may  be  regarded  as  responding  to  its  past,  present, 
or  future  experience  according  as  y < Q,  y = Q or  y > 0.  A different 
motivation  for  proving  that 


p > 0 

y < 0 


(2.12) 


can  be  found  in  Markovitz  5 Coleman  (1964). 

Although  the  special  solutions  of  the  second-order  fluid  are 
undoubtedly  useful  for  the  interpretation  of  the  behavior  of  real 
viscoelastic  fluids  under  appropriate  circumstances,  Coleman,  Duffin  5 
Mizel  (1965)  have  found  grave  defects  in  the  theory.  They  have  in  fact 
showed  that  all  solutions  for  unsteady  parallel  flow  of  a second-order 
fluid  are  unstable.  However,  it  is  not  certain  that  these  results 
render  this  fluid  useless.  For  example,  it  could  be  argued  that  as  the 
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disturbances  grow  unbounded  in  time,  they  are  excluded  from  the  class 
of  "slow"  flows  resulting  by  retardation  of  a given  flow. 


2.6  Oldroyd  Fluid 

Oldroyd  (1950)  in  a paper  entitled  "On  the  Formulation  of 
Rheological  Equations  of  State"  laid  down  what  is  now  known  as  the 
"Principle  of  Material  Indifference"  or  "Material  Objectivity":  "The 

form  of  the  completely  general  equations  must  be  restricted  by  the 
requirement  that  the  equations  describe  properties  independent  of  the 
frame  of  reference  . . . Moreover,  only  those  tensor  quantities  need  be 
considered  which  have  a significance  for  the  material  element 
independent  of  its  motion  as  a whole  in  space."  Oldroyd  in  effect 
observed  that  the  formulation  of  all  equations  in  terms  of  coordinates 
moving  and  rotating  with  the  local  vorticity  suffices  to  satisfy  this 
requirement . 

Let  coordinate  surfaces  ^ = constant  be  chosen  as  surfaces  which 
move  and  rotate  with  the  fluid  particles;  whether  in  addition  these 

J 

surfaces  also  deform  with  the  material  (i.e.,  are  embedded  in  the  fluid) 

will  not  be  specified  now.  In  these  convected  coordinates,  ^ and  time 

t are  taken  as  independent  variables;  for  example  a second  order  tensor 

k 

may  be  specified  as  6^(C,t).  Having  found  the  constitutive  equation  in 
terms  of  tensors  referred  to  their  convected  components  (i.e.,  the 
components  referred  to  a convected  system  of  reference) , for  the  sake 
of  simplicity  of  applications  it  is  necessary  to  express  all  the  tensor 
quantities  referred  to  an  arbitrary  curvilinear  fixed  coordinate  system 


1 
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x1,  having  a metric  tensor  of  g^^Cx) . Let  the  fixed  components  of 
k k 

6^(C,t)  be  b^x.t).  Now  in  the  constitutive  equation  written  in  con- 

v 

vective  frame,  there  may  have  been  partial  time  derivatives  3B^(£,t)/3t, 
holding  convected  coordinates  constant.  The  question  is,  how  do  these 
"convective  derivatives"  transform?  The  answer  is  provided  by  the 
following  theorem  of  Oldroyd: 

l 

'The  tensor  whose  convected  components  are  36"  *i' ‘ ‘ (C,t)/3t  has 
fixed  components 


k k 

d b".  3b"  . , , . 

c ..i..  ..i..  , ,nv..k..  „ m,..k..  „ ' k,  . .m. . 

yz = +UD.  + I u.D  +Ewb. 

dt  at  i..m..  m..i.. 


♦ Z eV'*" 

l . .m. . 


r 1 k,  . .m. . 
Z e b 

m . .1 . . 


(2.13) 


where  Z(Z  ) denotes  a sum  of  all  similar  terms,  one  for  each  covariant 

(contravariant)  suffix,  a comma  preceding  a suffix  denotes  covariant 

differentiation,  and  u1,  e^  denote  the  velocity  vector,  the 

corresponding  vorticity  and  rate  of  strain  tensors  (all  functions  of 

x1  and  t)  associated  with  the  motion  in  space  of  an  arbitrary  coordinate 

system  ^ , not  necessarily  of  the  fluid  unless  the  ^ system  is 

embedded  in  the  fluid."  Of  course  u>.,  = %(u,  •-  u.  ,)  and  e.,  = 

IK  K j 1 1 * K 1 K 

ui,k>- 

The  presence  of  dc/dt  in  a constitutive  equation  referred  to  fixed 
coordinates  obviously  conforms  with  the  principle  that  the  equation 
must  be  independent  of  arbitrary  time  dependent  rigid  motion.  Whether 
in  addition  the  coordinates  with  respect  to  which  the  original 
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constitutive  equation  was  written,  are  embedded  in  the  fluid  and  deform 

with  it  cannot  be  settled  by  recourse  to  invariance. 

In  case  the  t,  are  embedded  in  the  fluid,  Oldroyd  calls  d£/dt  the 

"convected  derivative"  and  denotes  it  by  6 / 5 1 . Then  u.  , e.  . and  ui.  . 

i ij  ij 

of  the  reference  frame  are  equal  to  the  corresponding  fluid  variables, 
i.e. , u.  = U.  , e.  . = E.  . and  w.  . = 12.  . . From  (2. 13)  therefore  the 

l l ij  ij  ij  ij 

convective  derivatives  of  second-order  covariant  and  contravariant 
tensors  are  of  the  form 


6b. . (x,t)  3b. . 

-rpi  ' = + UV  . + Um.b  . + Um.b. 

6t  3t  ij ,m  ,i  mj  , j lm 


6b^(x,t)  __  3b^  + yH'jjij  _ ui  bmj  . bim 

6t  ,m  ^in  ,m 


(2.14) 


If  the  frame  £ is  rigid  so  that  it  simply  moves  and  rotates  with 
the  fluid  but  does  not  deform  with  it  (e^  = 0),  then  Oldroyd  denotes 
d£/dt  by  V/Vt  and  calls  it  the  "material  derivative."  From  (2.15)  the 
material  derivative  of  b^  in  cartesian  tensor  notation  is 


Pb . . (x,t)  3b . . 3b . . 

- = — 5^-  + U,  — — + 0.,b,  . + 12..  b., 

Vt  3t  k 3x^  lk  kj  jk  lk 


(2.15) 


where  of  course  = '2(3U^/3x^  - Slh/Sx^). 

As  a matter  of  interest  it  may  be  mentioned  that  if  we  had  taken 
the  £ to  be  simply  translating  with  the  fluid,  neither  rotating  = 0) 
nor  deforming  (e^.=  0),  then  from  (2.13)  dc/dt  would  have  reduced  to  the 
familiar  D/Dt  which  Oldroyd  calls  the  "intrinsic  derivative": 
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Dbj(x.t) 

Dt" 


uV 

l,m 


To  summarize  in  rough  language,  D/Dt  is  a measure  of  the  rate  of  change 
which  allows  for  the  translational  motion  only,  V/Vt  makes  allowances 
for  the  rotational  as  well  as  the  translational  motion  of  the  element, 
and  6/fit  is  the  rate  of  change  which  takes  into  account  the  transla- 
tional, rotational  as  well  as  the  deforming  motion  of  the  fluid.  For  a 
scalar  these  three  derivatives  are  all  equal.  This  is  because  the 
value  of  a scalar  does  not  depend  on  any  coordinate  system,  and  hence 
D/Dt  of  pressure  or  temperature  denotes  an  intrinsic  property  of  a 
material  element.  However,  for  quantities  which  are  not  scalars  but 
are  measured  in  relation  to  certain  directions,  it  is  evident  that  to 
get  an  intrinsic  property  of  a material  -element  the  rate  of  change  must 
be  measured  in  some  way  related  to  directions  which  rotate  with  the 
material  element. 

Since  the  metric  tensor  g^  vanishes  when  operated  on  by  a/ at, 

D/Dt  or  V/Vt,  it  follows  that  the  operations  of  raising  and  lowering 

suffixes  .commute  with  these  types  of  differentiation.  On  the  other 

i k i k 

hand  it  can  be  shown  that  Sg^/fit  = 2E^,  fig  /fit  = - 2E  and  therefore 
the  operations  of  raising  or  lowering  a suffix  and  of  convected 
differentiation  do  not  commute  [if  they  did,  the  two  equations  in  [2.14) 
would  be  the  same] . 

Having  discussed  the  various  types  of  convective  derivatives, 
Oldroyd  illustrates  the  process  of  formulating  equations  of  state.  It 
has  been  found  that  the  behavior  of  certain  viscoelastic  liquids  at 
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small  rates  of  strain  can  be  approximately  described  by  a linear 
equation 


(1  ♦ A -£-)S.  . 

1 ij 


2u(l  + A„  £■) E.  . 

2 9t  lj 


(2.16) 


where  S. . is  the  extra  stress,  E. . = ^(U.  . + U.  .)  is  the  strain  rate, 
iJ  ij  i,J  J.i 


and  y is  the  viscosity;  A and  A are  called  the  "relaxation  time"  and 

1 2 


"retardation  time"  respectively,  with  the  direct  physical  significance 
that  if  the  motion  is  suddenly  stopped  the  stresses  decay  as  exp(-t/A^) 


and  if  the  stresses  are  removed  the  rates  of  strain  decay  as  expf-t/A^) 


Equation  (2.16)  was  shown  to  be  associated  with  a suspension  of  elastic 
spherical  particles  in  a Newtonian  viscous  liquid  by  Frohlick  § Sack 


(1946).  They  showed  that  the  introduction  of  3S^/9t  requires  a 


further  condition  because  at  certain  times,  e.g.,  on  application  or 


removal  of  stress,  S.,  will  be  allowed  to  be  a discontinuous  function 

ij 


of  time;  the  additional  condition  is  that  A S..-2yA  E..  is  continuous 

1 1J  2 1J 


at  a discontinuous  change  in  stress.  They  also  showed  that  the 
physical  model  is  such  that  A^  > A . It  is  apparent  that  as  A and  A2 
tend  to. equality,  the  material  must  become  more  and  more  exactly  a 
Newtonian  liquid  of  viscosity  y. 

In  addition  to  being  valid  only  at  small  rates  of  strain,  equation 
(2.16)  was  written  in  orthogonal  cartesian  coordinates.  It  is  thus  not 
for  general  use,  and  the  immediate  problem  is  that  of  generalizing  the 
equation  (i.e.,  formulate  an  invariant  equation)  so  that  the  equations 
of  state  can  be  applied  to  material  in  general  motion,  in  which  the 
rates  of  strain  will  not  be  necessarily  small.  Two  possible  generaliza- 
tions of  (2.16)  are 
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(1  ♦ X -£)S.  . = 2u(l  ♦ X ~)E.  • 

1 6t  lj  2 Ot'  1J 


(Oldroyd  Liquid  A) 


(2.17) 


where  S. . and  £.  . are  covariant  absolute  tensors,  and 

ij  ij 

(1  + A^  ^)S1J’  = 2u(l  ♦ A2  ^-)E1;i  (Oldroyd  Liquid  B)  (2.18) 

where  S1-^  and  E1-1  are  contravariant  components.  Lumley  (1971a)  showed 
that  if  the  polymer  molecules  are  represented  by  dumbell  models,  then  a 
dilute  solution  of  polymer  molecules  obeys  (2.18)  with  A = T,  A^  = 

T/(l  + c[m]),  u = n ( 1 + c[y]),  where  T is  the  terminal  relaxation  time 
of  the  molecules  (that  is,  the  initial  small  departure  of  the  distribu- 
tion from  spherical  symmetry  relaxes  to  zero  with  time  constant  T) , [u] 
is  the  intrinsic  viscosity,  c is  the  concentration  of  polymer  and  n is 
the  viscosity  of  the  solvent. 

Oldroyd  (1958,  1961)  showed  that  the  most  complete  generalization 
of  (2.16)  is,  in  cartesian  coordinates. 


VS 


S.,  + A 


ik 


+ p S..E.,~u  (S , . E . , +S  . , E . . ) + v S./>E.p5.. 


ik  T q j j ik  i ^ i j j k j k i j J i ik 


= 2w 


PE.. 

n , 

E . , + A — TyT 

lk  2 Pt 


2y  E . .E  + v E . »E . »6  . , 
2 iJ  jk  2 ik 


(2.19) 


where  u , p , p , v and  v are  five  more  arbitrary  constants,  each  with 
0 12  1 2 

the  dimensions  of  time.  The  generalization  (2.19)  contains  terms 
linear  in  the  stresses  and  quadratic  in  stresses  and  space  derivatives 
of  velocity  taken  together;  other  generalizations  of  (2.16)  must  contain 
terms  of  higher  degree  in  the  stresses  and  velocity  gradients.  The 
Oldroyd  liquid  A,  represented  by  (2.17),  is  a special  case  of  (2.19)  with 
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X = y > X = y >0 
112  2 


y = V = v =0 
0 1 2 


It  has  been  shown  by  Oldroyd  (1950)  that  liquid  B would  exhibit  the 
positive  Keissenberg  effect,  whereas  liquid  A would  exhibit  the  negative 
Weissenberg  effect  (i.e.,  sinking  down  near  the  inner  cylinder). 

Oldroyd  (1958,  1961)  has  shown  that  the  constants  in  (2.19)  have 
the  following  significance  [quoting  from  Oldroyd  (1961)]: 

"(a)  The  behavior  at  small  rates  of  strain  is  realistic,  provided 

y > 0 

X > X >0 
1 2 

"(b)  The  variable  apparent  viscosity  is  a decreasing  function  of 

the  rate  of  shear,  varying  from  the  limiting  viscosity  y at  low  rates 

of  shear  to  a lower  limit  a y/a  at  high  rates  of  shear,  provided 

2 1 r 


a>a>7ro>0 

1 2 ~ 9 1 


where 


a = X + y (y  - t v ) ■ y (y  - v] 
1 1 0 1 2 1 11  1 


a = XX  + y y - Tv)-y(y  - v) 
2 12  0 2 ^ 2 12  2 


"(c)  The  normal  stress  [effect  corresponds  to]  a simple  tension 
along  the  streamlines,  in  excess  of  the  isotropic  state  of  stress  in 
directions  normal  to  the  streamlines,  in  steady  shearing  flow  between 
a cone  and  a flat  plate,  provided 
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U = X 
1 1 


u = X 
2 2 


2 2 


"(d)  A liquid  whose  constants  satisfy  (a),  (b)  and  (c)  will 
necessarily  show  a positive  h'eissenberg  climbing  effect,  which  is  that 
most  often  observed." 

Combining  (a),  (b)  and  (c) , we  are  led  to  the  following  restric- 
tions on  the  constants: 

y > 0 

A > X > 0 
1 2 


Ul=  \ 


U = X 
2 2 


(2.20) 


yX+(A--flJ)u>uA+(A-|-y)\;  — 4[d  X + (x  - -|u)v]>0 

01  1 2 q 1 0 2 1 2 o 2 y01  1 2 o 1 


(X  - y ) (A  v - X v ) > 0 
1 2 o 1 2 2 1 


Observe  that  any  one  of  the  constants  y , v and  v may  vanish  con- 

0 12 

sistently  with  these  conditions,  but  if  two  vanish  the  only  interesting 
possibility  is  given  by 

y > 0 


y > 0 
0 


X=y>X=y>  — X >0 
1 1 2 2 “ y 1 


(2.21) 


v = v = 
1 2 


0 
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which  differ  from  the  liquid  B only  in  the  presence  of  a p term  in  the 

0 

equation  of  state. 

A few  comments  must  be  made  about  the  experimental  measurements  of 
the  elastic  constants  in  the  Oldrovd  equation.  Toms  5 Strawbridge 
(1953)  measured  the  A ^ and  A ^ of  some  polymer  solutions  in  an  indirect 
way.  They  supposed  that  the  linearized  Oldroyd  equation  (2.16)  is  valid 
for  small  rates  of  strain,  and  obtained  the  theoretical  solutions  for 
the  amplitude  of  oscillation  of  an  inner  cylinder  (of  known  moment  of 
inertia  and  held  by  a torsion  wire  of  known  spring  constant)  when  the 
outer  cylinder  oscillates  at  a known  amplitude  and  frequency.  They 
plotted  these  solutions  as  a graph  of  9,  defined  as  the  ratio  of  the 
amplitudes  of  the  inner  and  outer  cylinders,  against  frequency  w,  with 
X^  and  regarded  as  parameters.  Then  they  obtained  the  experimental 
0-versus-oj  curve  of  a polymer  solution,  and  from  a comparison  of  the 
experimental  and  theoretical  plots  arrived  at  approximate  values  of  X 
and  X^  of  the  solution.  The  important  conclusion  they  arrived  at  is 

that  while  the  viscosity  u and  the  stress  relaxation  time  X increased 

1 

steadily  with  the  concentration  of  the  polymer,  the  strain  rate 
relaxation  time  A^  on  the  other  hand  remained  independent  of  the  con- 
centration. They  suggested  that  perhaps  different  molecular  mechanisms 
are  responsible  for  causing  A and  X in  a solution. 


2.7  Walters  Fluids 


Let  Ay^j(t  ) be  the  strain  at  a stationary  element  of  material  at 
time  t , which  causes  an  extra-stress  (t)  at  the  same  point  at  ti 
t.  Let  us  assume  that  the  cause  and  effect  are  related  by 


me 


As..(t)  = iKt-t')AYij  Ct ' ) 


where  tp  may  be  called  the  influence  function  or  relaxation  function. 
The  above  equation  simply  states  that  the  material  under  consideration 
possesses  a memory  for  its  past  history. 

Suppose  the  material  is  subjected  to  a series  of  step  changes  in 
the  strains.  The  resultant  stress  at  time  t can  be  written  as 


Sij(t)  = <Kt-t  3 AYij  Ct  ) + <|<(t-t  )Ayij  (t  ) + ... 

if  linear  (Boltzmann)  superposition  is  valid.  If  the  sequence  of  strains 
experienced  by  the  material  point  is  not  a series  of  step  functions  but 
rather  a continuously  varying  function  of  time,  we  replace  the  summation 


by  an  integration  to  obtain 


f , dY  , , 

S±  - (t)  = Ht-t  ) — il  (t  )dt 

J dt 

-00 

» » » 

Denoting  dy„  (t  ) /dt  = 2E^  (t  ),  the  strain  rate  tensor,  we  get 


Si:j(x,t)  = 2 <Kt-t  3E.j(x,t  ) dt 


as  the  extra-stress  at  a stationary  element  of  material  at  x at  time  t. 
Walters  (1960)  has  arrived  at  the  above  result  by  superposing  an 
infinite  number  of  Maxwell  elements  in  parallel.  If  N(t)dT  represents 
the  total  viscosity  of  all  the  Maxwell  elements  with  relaxation  times 


between  t and  t + dt,  then  he  derives  (2.22)  with  an  influence  function 


e dx 
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N(t)  is  called  the  distribution  of  relaxation  times,  or  the  relaxation 
spectrum,  a characteristic  of  the  material. 

Quoting  from  Walters  (1962a),  note  that  in  a general  flow  " ...  the 
equation  of  state  (2.22)  must  be  associated,  not  with  the  particular 
point  x*  in  space  over  a period  of  time,  but  with  the  same  material 
element  over  a period  of  time:  the  element  which  is  instantaneously  at 
the  point  x1  at  time  t.  In  consequence,  the  integral  occurring  in  (2.22) 
must  bo  regarded  as  a convected  integral  following  the  moving  material. 
Only  under  the  restrictive  conditions  of  small  rates  of  shear  in  a sta- 
tionary material  element  can  the  convected  integral  be  written  as  in 


(2.22). 

"The  problem  is  now  one  of  generalising  equation  (2.22)  into  a 
form  that  can  be  applied  to  materials  in  general  motion,  in  which  the 
rate  of  strain  is  not  necessarily  small.  This  of  course  involves 
replacing  the  integral  occurring  in  (2.22)  by  a convected  integral 
following  the  moving  material.  Oldroyd  (1950)  has  included  transforma- 
tion of  convected  integrals  in  his  analysis."  He  proved  that  any 
t i » 

integral  t f 3 (S,t  )dt  transforms  into 

tQ  . . j . . 


i i i i i 

when  referred  to  fixed  coordinates,  wrhere  x = x (x,t,t  ) is  the 

position  at  time  t'  of  the  element  which  is  instantaneously  at  the  point 

* £ 
x1  at  time  t,  and  3 ’ ' ‘ ' are  the  convected  components  of  b*\’*. 

• • J • • • • j • • 

Walters  (1962a)  points  out  that  "in  general,  the  displacement  functions 

I i 

x cannot  be  expressed  as  simple  functions  of  the  usual  kinematic  vari- 
ables and  it  is  necessary  to  determine  them  separately  in  each 
individual  flow  problem." 
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Using  Oldroyd's  rule  and  noting  that  the  suffixes  may  be  raised  in 
(2.22)  under  the  restrictive  conditions  in  which  the  equation  is 
supposed  true,  the  two  simple  possible  generalizations  of  (2.22)  are 


- * in  1 r 

Si^(x,t)  =2  / tKt-t  ) — - — — j—  Emr(x  )dt  (Walters  Liquid  A ) 

(2.24) 

t ■»  V 

t I 2*  JW  TTrr  f * 1 * 

2 / i|/(t-t  ) — i — — — i — E (x  ,t  )dt  (Walters  Liquid  B ) 

„ ra  r 

(2.25) 


Slk(x,t) 


3x  ax 


■ I 

The  designations  "Liquid  A " and  "Liquid  B " were  given  by  Walters 
(1962a),  and  have  been  used  by  other  writers.  The  liquids  designated 

A and  B by  Oldroyd  (see  equations  2.17  and  2.18)  are  special  cases  of 

• » 

A and  B respectively,  obtained  by  substituting 

\ A -A 

N(t)  = v 6(x)  ♦ p -L-±  6(t-X  ) 

i 1 * 


in  equations  (2.23),  (2.24)  and  (2.25).  The  Maxwell  liquid  with  one 
relaxation  time  at  t = is  also  a special  case  obtained  by  writing 
N(t)  = p6(r-A1),  while  the  Newtonian  liquid  is  obtained  by  taking  N(t)  = 

p6(i). 

t f 

Some  of  the  simplest  consequences  of  equations  of  state  A and  B 
have  been  investigated  by  Walters  (1962a)  for  a number  of  steady  flow 
problems  involving  simple  shearing.  No  difficulty  is  encountered  in 
such  problems  since  the  displacement  functions  x'1  can  be  obtained 
without  difficulty  and  the  rate  of  strain  E^  or  E in  any  element 

I 

does  not  vary  with  t and  so  can  be  taken  outside  the  integrals. 
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Walters  (1962b)  also  considered  simplifications  of  his  liquids  A 
and  b'  for  the  special  case  of  short  memory.  The  derivation  of  his 
simplifications  is  carried  out  by  referring  the  equations  to  a convected 
coordinate  system  ^ drawn  in  the  material  and  deforming  continuously 
with  it.  In  the  case  of  liquid  A*,  the  equation  of  state  obviously 
becomes 


o-jCg.t)  = 2 / <Kt-t  )eij(§,t  )dt 


where  a. . and  e. . are  the  convected  components  of  S. . and  E. . respec- 
ij  ij  ‘ ij  ij 

tively.  If  the  rate  of  strain  tensor  is  now  expanded  in  a Taylor  series 
about  the  time  t,  we  have 

Eij(§,t')  = e^fg.t)  - (t-t')  ^-e.jCg.t)  + ... 


where  3/3t  denotes  a partial  differentiation  with  respect  to  time 
holding  convected  coordinates  constant.  The  constitutive  equation 
then  reduces  to 

t 3e.  * t 

a = 2e . . / (Kt-t')dt'  - 2 / (t-t ' ) i|»(t-t ' ) dt ' ♦ 


Using  (2.23)  and  neglecting  terms  involving  Qf  tnN(t)dT,  (n  >.  2) , this 
can  be  written  as 


o . . = 2 pc . . - 2K  e . . 

ij  ij  3t  ij 

00 

where  u = / N(t)dT  = limiting  viscosity  at  small  rates  of  shear,  and 

00 

K = / TN(x)dt.  Transforming  the  constitutive  equation  from  convective 

to  a fixed  coordinate  system  x1,  we  finally  obtain 
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2pE. . - 2K 
ij 


(Walters  Liquid  A 
with  short  memory) 


(2.26) 


where  6/5t  is  Oldroyd's  convective  derivative. 

In  a similar  manner  it  can  be  shown  that  the  simplified  form  of  the 
equation  of  state  for  liquid  B in  case  of  short  memory  is 


S1^  = 2uE1J'  - 2K  ^1— 

ot 


(Walters  Liquid  B 
with  short  memory) 


(2.27) 


It  is  interesting  to  note  that  we  have  found  that  equations  (2.26)  and 
(2.27)  are  also  derivable  from  Oldroyd’s  equations  for  slow  flows.  Let 
us  illustrate  for  the  Oldroyd  liquid  A: 


6S.  . 

S.  . + A —ril 
1J  I St 


6E.  . 

= 2u(E..  + A -A 

1J  2 


Now  for  slow  flows  (or  small  viscoelasticity,  or  short  memory)  the 
terms  multiplied  by  A and  A will  be  smaller  compared  to  the  other  two 
(Newtonian)  terms.  Thus,  no  great  error  will  be  introduced  if  in  the 
term  multiplied  by  A^  we  substitute  S^_.  = 2uE^ ^ . Then  Oldroyd's 
liquid  A becomes 

6E . . 

S. . = 2yE. . - 2p(A  -A  ) 

ij  ij  1 2 6t 


which  is  identical  to  (2.26). 


2.8  Normal  Stress  Functions  in  Viscoelastic  Fluids 

Coleman  5 Noll  (1961)  have  shown  that  a general  memory  fluid 

• * 

( simple  fluid  ) displays  the  following  stress  pattern  in  a simple 
shearing  flow  of  shear  rate  F: 


S12  = n(r)r 

S-S  = o (r)  S = S =0  (2.28) 

11  33  1 13  23 

S-S  = a (H 
22  33  2 

where  the  direction  1 is  taken  along  the  flow  and  direction  2 is  taken 
perpendicular  to  the  plates.  n(T)  is  the  shear  dependent  viscosity,  and 
a.  and  a.  are  the  so-called  first  and  second  normal  stress  differences, 


respectively.  The  three  material  functions  n(F),  0^(0  and  a (T)  are 
all  even  functions  of  the  shear  rate  T because  of  isotropy,  and  of 
course  0^(0)  = a2(0)  = 0 from  the  hydrostatic  condition;  the  three 
functions  are  otherwise  quite  arbitrary. 

Another  definition  of  the  first  normal  stress  difference  quite 
widely  used  [see,  for  example,  Mclntire  5 Schowalter  (1970),  Ginn  5 


Denn  (1969)]  is 


N (r)  = S - S 

1 11  22 


(2.29) 


There  have  been  a number  of  attempts  to  measure  the  normal  stress 
functions  of  polymer  solutions.  The  experimental  difficulties  are, 
however/  severe  and  there  is  considerable  disagreement  in  the  results. 
Markovitz  (1963)  found  that  the  Reiner-Rivlin  result  o = o2  and 
the  Weissenberg  conjecture  = 0 are  both  in  variance  with  their 
experimental  data;  he  also  found  that  > °2  > 0 in  many  cases. 

Denn  6 Roisman  (1969)  found  that  N is  always  positive  and  is  always 
negative.  A negative  was  also  reported  by  Ginn  (1968).  however, 
Huppler  (1965)  reported  a which  are  not  only  positive  but  also  far 


52 


greater  in  magnitude,  by  factors  of  ten  or  more,  than  the  results  of 
Denn  5 Roisman  (1969). 

The  theoretical  predictions  of  the  normal  stress  functions 
according  to  the  commonly  used  constitutive  equations  are  definitely  of 
interest.  Consequently  we  have  made  a list  of  them,  which  is  displayed 
in  Table  2.1. 


2.9  Normal  Stresses  and  Weissenberg  Effect 


During  working  with  the  commonly  used  constitutive  equations,  we 
came  to  the  observation  that  the  sign,  and  magnitude,  of  the  Weissenberg 
effect  may  depend  on  the  method  of  observation,  namely  through  the  shape 
of  the  free  surface  or  through  the  height  difference  of  radial  tappings. 
This  will  be  explained  below. 

(i)  An  Apparent  Contradiction 

Let  us  consider  three  well-known  constitutive  equations,  the 
Oldroyd  A fluid 


6S . 


S.  . + X 

ij  1 


11 


5t 


2“<Eij 


X > X > 0 
1 2 


(2.30) 


the  Walters  A*  fluid  with  short  memory 

6E.  . 

S.  . = 2uE . . - 2K  — s-ii  , K > 0 (2.31) 

lj  ij  6t 


and  the  second-order  fluid 


5A . . 

S.  . = pA.  . + SA.A.  . + y ■ « 

lj  lj  l kj  6 1 


Y < 0 


(2.32) 
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TABLE  2.1 

Viscometric  Functions  of  Common  Fluids 


Fluid 


Second-Order 


Walters  A 


Equation 


Normal  Stresses  in 

Viscometric 

Simple  Shear 

Functions 

sn*  6rZ  2 

o = er2 

S22=  (e+2y)r 

o » (B+2y) r 

S * 0 

N = -2Yr 

33 

1 

S = 0 

a = 0 

11 

1 

S = -2T  r xN(x)dx 
22  0 ^ J 

a2=  -2r  Q/  xN(x)dT 

2 00 

s = o 

33 

N = 2rQ/  xN(x)dx 

Walters  B 


Eq.  (2.25) 


Walters  A 
(short  memory) 

Eq.  (2.26) 

S = 0 

Li  2 

S = -2Kr 
22 

S 3 3=  ° 

a = 0 

a = -2Kr2 
nx=  2Kr2 

Walters  B 
(short  memory) 

Eq.  (2.27) 

sir  2Kr2 
S22=  0 
S33=  0 

a^=  2Kr2 

a = 0 

L 2 
Nj=  2Kr 

Oldroyd  A 

Eq.  (2.17) 

sn=  0 

s22=  -2u(-\1-X2)r2 

S = 0 

3 3 

V 0 

a2=  -2u(X1-X2)r 
N j=  2u(Xx-X2)r2 

Oldroyd  B 

Eq.  (2.18) 

sn=  2u(xrx2)r2 

S22-  0 
S33=  0 

ar  2y cx1~x2)r2 

V 0 2 

V 2^YVr2 

Fluid 

Constitutive 

Equation 

Normal  Stresses  in 
Simple  Shear 

Viscometric 

Functions 

Maxwell  A 

Eq.  (2.58) 

sn=  o , 

S = -2pA  r 
22  1 

S = 0 

33 

V 0 2 

02=  - 2u A ^ r 
N1=  2uAir2 

s = 2ux  r2 

2 

a = 2uX.r 

ii  i 

1 1 

Maxwell  B 

Eq.  (2.59) 

S22=  ° 

a = 0 

* 2 

S = 0 

N = 2pX,r 

33 

1 1 

2 

2 

s = f r 

a = for 

11  2 2 

1 * 2 

Reiner-Rivlin 

Eq.  (2.7) 

S22=  f2r 

a2-  f2r 

S = 0 

N = 0 

33 

1 
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where  we  are  using  general  tensor  notation  and  the  5/ 5 1 in  each  of  the 
three  equations  is  given  by  the  first  equation  of  Eq.  (2.14).  It  is 
apparent  that  a result  for  the  Walters  A1  fluid  (with  short  memory) 
could  be  deduced  either  by  substitu..  ig  X^  = 0,  uX2  = -K  into  the 
corresponding  result  for  Oldroyd  A fluid,  or  by  substituting  S = 0, 
y = -K  into  a result  for  the  second-order  fluid.  Using  the  existing 
results  for  the  Oldroyd  A and  second-order  fluids,  let  us  therefore 
predict  the  sense  of  the  Weissenberg  effect  for  the  Walters  a'  fluid 
(with  short  memory)  in  circular  Couette  flow  - positive  (climbing  near 

inner  cylinder)  or  negative  (dipping  down  near  inner  cylinder). 

12  3 

Let  x , x and  x denote  the  cylindrical  coordinates  r,  6 and  z 
respectively,  with  z taken  along  the  axis  of  the  cylinders.  The 
physical  components  of  all  tensors  will  be  denoted  by  using  parentheses 
around  the  suffixes.  Markovitz  5 Coleman  (1964)  have  shown  that 
ignoring  inertia,  the  difference  of  physical  components  of  compressive 
radial  stresses  between  two  points  at  the  same  horizontal  level 

z but  at  radial  positions  and  (>  Rp  is  given  by 


(ID 


YC 


(2.33) 


where  terms  in  3 dropped  out  and  C is  a positive  constant.  Since  y < 0, 
they  conclude  that  the  second-order  fluid  displays  the  positive  Weissen- 
berg  effect.  And  since  6 has  automatically  vanished,  it  also  follows 
that  the  same  conclusion  must  also  be  true  of  the  Walters  A'  fluid  with 
short  memory. 

On  the  other  hand  Oldroyd  (1950)  has  shown  that  for  the  Oldroyd  A 


fluid 
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R2  ^1  » 

-TC33)  + T(33)  = <W  C (2-34) 

I 

where  C is  a positive  constant.  Since  X^  > X , it  follows  that  the 
Oldroyd  A fluid  displays  the  negative  Weissenberg  effect.  (Since  a 
real  viscoelastic  fluid  generally  displays  the  positive  Weissenberg 
effect,  Oldroyd  therefore  discards  this  constitutive  equation  and 
accepts  the  Oldroyd  B fluid  which  predicts  a positive  Keissenberg 
effect.)  And  by  substituting  = 0,  X^  = -K/y  into  Eq.  (2.34),  it 
follows  that  the  same  conclusion  must  also  be  true  of  the  Walters  A* 
fluid  with  short  memory. 

This  reason  for  this  apparent  contradiction  regarding  the  sense  of 
the  Weissenberg  effect  of  the  Walters  A*  fluid  (with  short  memory)  will 
now  be  investigated. 

(ii)  Consideration  of  Oldroyd  A Fluid 
Let  the  physical  components  of  velocity  be 


U0)  ' 0 

U,  . = no 

(2) 


“(3,  ■ ° 


(2.3S) 


From  Eq.  (2.30)  it  can  be  deduced  that  the  physical  components  of  extra 
stress  are 


“W 


-2u(Xr»2)(I^)2 


dm  n 

Pr  d?  ° 


yr 


dio 

dr 


(2.36) 


0 


0 


0 
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The  equation  of  motion  is 


ci  ii  cij 

= f - P Hg  + S J 

> J » J 


(2.37) 


where  f1  is  the  body  force  and  g1^  is  the  metric  tensor.  Using  (2.36), 
an  integration  of  the  equation  of  motion  in  the  6-direction  gives  the 


velocity  distribution 


a)  = a 


(2.38) 


where  a and  b can  be  found  in  terms  of  the  angular  velocities  and  radii 
of  the  boundaries.  IVith  Eq.  (2.38),  an  integration  of  the  equation  of 
motion  in  the  radial  direction  yields  the  pressure  distribution 


P = P(h  a2r2 — - 2ablogr)  - 6p  (A  -A  )^— 

2r  1 2 r 


(2.39) 


The  effect  of  inertia,  represented  by  the  first  term  on  the  right  hand 
side  of  the  above  equation,  is,  of  course,  to  lower  the  free  surface 
near  the  inner  cylinder.  To  find  out  the  sense  of  the  Keissenberg 
effect,  one  has  to  ignore  this  term  and  concentrate  on  the  viscoelastic 
term.  The  compressive  stress  on  horizontal  planes  is  - T,  . = P - 
S(33)’  311  ^ since  = it  follows  from  Eq.  (2.39)  that,  ignoring 


inertia, 


R R b (R  - 

- T(33)  + T(33)  " 6^YV  "77 


2 4 4 

b (R  -R  ) 
2 r 


(2.40) 


which  shows  that  the  free  surface  will  dip  near  the  inner  cylinder. 

However,  let  the  observation  be  carried  out  through  radial  holes,  as 

shown  in  Fig.  2.1.  Using  -T.  ^ = P - Eq.  (2.39),  and 

2 4 

^(11)  = ^(A^-A^b  /r  > it  follows  that 


Eq.  (2.39),  and 
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R2  R1 
T «f  T 

(n)  do 


- 2p  (X  -X  ) 
v 1 2 


2 4 4 

b (R  -R  ) 

2 V 

_ 4 4 

R R 
1 2 


(2.41) 


which  shows  a positive  We issenberg  effect. 

The  explanation  of  f.e  apparent  contradiction  mentioned  before  is 

I 

now  clear:  The  Oldroyd  A fluid,  as  well  as  the  Walters  A fluid  with 

short  memory,  display  the  negative  Weissenberg  effect  as  far  as  the 
free  surface  is  concerned,  but  they  would  display  the  positive  Weissen- 
berg effect  when  observed  through  radial  tappings. 

As  a matter  of  interest,  we  may  note  that  for  the  Oldroyd  B fluid, 
which  is  widely  accepted  as  one  of  the  most  reliable  of  all  simple  con- 
stitutive equations,  it  is  straightforward  to  show  that  the  Weissenberg 
effect  is  positive,  and  equal,  in  both  types  of  measurements.  It  can 

t I 

also  be  shown  that  in  a circular  Couette  flow  the  Walters  A and  B 
fluids,  without  short  memory,  behave  exactly  like  the  Oldroyd  A and  B 
fluids  respectively. 

(iii)  Consideration  of  Second-Order  Fluid 

Although  the  apparent  contradiction  is  already  resolved,  a few 
more  points  may  be  noted  from  a consideration  of  the  second-order  fluid, 
for  which  it  can  be  easily  shown  that 


] 


(3+2y) (r 


do).  2 

dr"^ 


0 


ur 


dw 

dr 


6^2 


(2.42) 


0 


0 


0 
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2 2 

P = p(J*  a r - ~ 
2r 


2ablogr) 


2(26+3y)b 


- Pgz 


(2.43) 


(2.44) 


The  difference  of  radial  thrusts  is  then  found  to  be 


r2  R 

T + T 

(11)  (11) 


2 4 4 

2yb  (R2-Rt) 

r,4_  4 

R R 
1 2 


< 0 


(2.45) 


which  shows  a positive  IVeissenberg  effect.  But  if  the  effect  is 
defined  in  terms  of  the  free  surface,  then  at  a particular  depth  it  is 
found  that 


2 4 4 

R R b (R  -R  ) 

- T,  , + T,  , = - 2(23+3y)  — — 

(33)  (33)  ^ _ 4„ 4 


(2.46) 


R R 
1 2 


which  can  be  of  either  sign,  depending  on  the  relative  values  of  3 and 
Y and  would  represent  a positive  IVeissenberg  effect  if  23  + 3y  > 0 and 
a negative  IVeissenberg  effect  if  23  + 3y  < 0. 

Markovitz  (1963)  has  measured  y of  a fluid  using  Eq.  (2.45),  to 
which  of  course  a correction  due  to  the  presence  of  inertia  was  made. 
They  then  measured  (3S+4y)  of  the  fluid  from  the  normal  stress  measure- 
ments in  steady  torsional  flow  between  parallel  plates,  and  thus  arrived 
at  the  values  of  3 and  y separately.  However,  it  appears  to  us  that, 
instead  of  doing  the  torsional  flow  experiment,  the  (23+3y)  of  the 
fluid  can  be  easily  measured  from  Eq.  (2.46).  Thus,  both  3 and  y of 
the  fluid  can  be  simply  derived  from  the  measurements  of  AT^^  and 
AT^3)  of  Fig.  2.1.  Of  course  the  normal  stress  measurements  in 
torsional  flow  are  no  more,  and  no  less,  reliable  than  the  measurement 
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of  the  free  surface  heights  of  the  circular  Couctte  flow.  However, 
there  is  an  important  point  in  favor  of  the  method  we  are  suggesting. 
This  is  the  fact  that  the  velocity  distribution  of  a simple  torsional 
flow  satisfies  the  equations  of  motion  only  if  inertia  is  neglected, 
whereas  the  circular  Couette  flow  satisfies  the  equations  of  motion 
irrespective  of  the  inertia,  and  a correction  for  inertia  can  always 
be  made  to  Eq.  (2.46). 

Alternatively,  if  the  fluid  is  not  assumed  to  follow  any  particular 

constitutive  equation,  then  the  experiment  we  are  suggesting  may  be  used 

to  determine  the  normal  stress  functions  a (T)  = S,  . - S,  . and 

1 (11)  (33) 


o2cn 


(22) 


of  a general  memory  fluid.  There  is  considerable 


controversy  regarding  even  the  sign  of  these  quantities  (see  Section 

2.8) . 


(iv)  Conclusions 

That  a constitutive  equation,  or  a real  fluid,  displays  the 
positive  or  negative  Weissenberg  effect  is  not  a sufficient  statement. 
The  method  of  observing  the  effect  - through  the  shape  of  the  free 
surface  or  through  the  height  difference  of  radial  tappings  - must  be 
mentioned.  The  Oldroyd  A fluid  and  the  Walters  A*  fluid  display  the 
negative  Weissenberg  effect  as  far  as  the  free  surface  is  concerned, 
but  they  would  display  the  opposite  effect  when  observed  through  radial 
tappings.  The  Oldroyd  B fluid  as  well  as  the  Walters  b'  fluid  display 
identical  positive  effect  in  both  types  of  observations.  The  second- 
order  fluid  would  display  the  positive  Weissenberg  effect  when  observed 
through  radial  tappings,  but  it  may  display  either  effect  when  the  free 
surface  is  observed.  A simple  experiment,  namely  that  of  observing  the 


Ik 
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Weissenberg  effect  by  the  two  methods,  is  suggested  for  measuring  the 
normal  stress  functions  of  a viscoelastic  liquid. 

2.10  El  astic  Energy  of  Viscoelastic  Liquids 

For  a purely  viscous  liquid  all  the  work  spent  in  deforming  a 
material  particle  is  dissipated  into  internal  thermal  energy.  It  is 
expected  that  in  a viscoelastic  liquid  a part  of  the  shearing  deforma- 
tion work  must  be  stored  as  elastic  potential  energy.  That  this  must 
be  the  case  is  evident  from  the  recoil  exhibited  by  such  liquids  upon 
release  of  stress.  Although  quite  a large  amount  of  literature  exists 
on  the  general  thermodynamics  of  continua,  no  work  seems  to  have  been 
done  on  the  specific  problem  of  finding  expressions  for  the  dissipation 
and  elastic  storage  of  the  various  constitutive  equations.  In  the 
remaining  sections  of  this  chapter  we  shall  make  an  attempt  in  this 
direction . 

The  differential  form  of  the  first  law  of  thermodynamics  for  any 
fluid  motion  is  [for  derivation  see  for  example  the  book  of  Aris  (1962) , 

p.  121] 


De 


p = T.  .U.  . - q.  . + 
Dt  lj  i,j  Hi,x 


pr 


(2.47) 


where  p is  the  density,  T„  is  the  stress  tensor,  q^  is  the  conduction 

heat  flux  vector,  r is  the  rate  of  heat  addition  (e.g.,  by  radiation) 

per  unit  mass,  and  e is  the  internal  energy  per  unit  mass;  for  our 

purposes  e should  include  both  thermal  and  elastic  energies.  T. .U.  . is 

b ij  i,J 

the  work  done  by  stresses  in  deforming  the  fluid  particles  per  unit 
volume  per  unit  time;  we  shall  call  it  the  deformation  work . For  an 
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incompressible  Newtonian  fluid  the  term  T„IL  ^ can  be  easily  shown  to 
be  equal  to  the  product  of  the  viscosity  and  a positive  quantity 
depending  on  the  strain  rates,  so  that  the  entire  amount  of  deformation 
work  goes  to  viscous  dissipation.  For  a viscoelastic  fluid  the  rate  of 
deformation  work  is  expected  to  be  the  sum  of  the  rate  of  dissipation 
and  the  rate  of  increase  of  elastic  strain  energy. 


2.11  Elastic  Energy  of  Second-Order  Fluid 


The  constitutive  equation  of  a second-order  fluid  is 


DA.  . 

Tij  ’ * “Aij  * 6AikAkj  * * Aikuk,j  * Vk.i’ 


where  we  have  dropped  the  subscript  1 on  the  first  Rivlin-Ericksen 

tensor  A. . = U.  . ♦ U.  ..  Multiplying  U.  . and  summing,  the  deformatic 
ij  i.J  * } 5 i,J 


work  is 


T.  .U.  . = -PU.  . + uA. .U.  . + BA..  A,  .U. 
ij  iJ  i.i  ij  i,J  ik  kj  i,j 


+ Y 3t  Ui,j  + UkAij,kUi,j  + (AikUk,j  + AjkUk,i^Ui,j 


(2.48) 


Now  let  us  simplify  the  various  terms  of  the  above  equation.  The  second 
and  third  invariants  of  A. . will  be  denoted  by  II  = A. .A. . and  III  = 

ij  ij  ij 

A^Aj^A^  respectively. 


A. .U.  . = h A. .A. . = h II 


ij 


ij  iJ 


klk\iUi,i  - AikAkj  «UiJ  + Uj,i)  = % AikAkjAij  = * 111 
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3A. . . _TT 

- i - u ^ (-a  \ a _ j,  3 1 1 

" 2*.  A..  " '•*  TT  tA.  - A.  .)  - <j  7— 


3t  ij  “ 3t  ij  •*  3t  v ij  ij- 


ukAij,kuij  - % ukAijfkAij = * uk<AijV.k = * ukn,k 


(AikUk,j+AjkUk,i)UiJ=^AikUkJ+AjkUk,PAij=  ^ ^AikAi j ^k , j +Aj  kAi j ^k , i ^ 


= h (A.,  A.  .A,  . + A. .A.  .A..)  = % III 

v ik  ij  kj  jk  ij  ki' 


With  these  simplications , (2.48)  reduces  to 


Vi.j-  “ ^ * 6 i|i  . ,0  ||i  * i ukn_t  . % m) 


which  can  be  rewritten  as 


T.  .U.  . = Ys  + £ 511 
ij  i,j  S 4 Dt 


(2.49) 


where 


fs  = j II  + %(B+y)III 


(2.50) 


Equation  (2.49)  has  a striking  similarity  with  the  corresponding 
equation  for  a Newtonian  gas  which  can  store  deformation  work  as 
"elastic  energy"  of  compression.  (Note  that  for  an  incompressible 
second-order  fluid  we  are  talking  of  elasticity  of  shape , whereas  for  a 
Newtonian  gas  we  are  talking  of  elasticity  of  volume.)  For  a Newtonian 


T.  . = - P5. . + pA. . - £ A. . 5. . 

ij  ij  3 kk  ij 


it  can  be  easily  shown  that 


T-  -U-  • 
lj  i.J 

Rate  of 

Deformation 

Work 


N 

Rate  of 
Dissipation 


P Dp_ 
p Dt 

Rate  of  Increase 
of  "Elastic" 
Energy 


(2.51) 


¥ s -d.  II  - — I 

N 2 1 6 

and  I = A. . is  the  first  invariant  of  A. ..  It  can  be  easily  shown  that 

n ij 

the  above  can  be  rewritten  as 


2 2 2 2 2 2 
= p(A  +A  +A  ) + rt(A  -A  ) + (A  -A  ) + (A  -A  ) ] 

N 12  23  31  6 11  22  22  33  33  11 


so  that  is  non-negative.  Equation  (2.51)  has  the  straightforward 
interpretation  written  below  it;  the  elastic  energy  of  compression  would 
result  from  a density  rising  cither  with  time  or  downstream,  as  in  the 
steady  flow  through  a diffuser.  Comparing  (2.49)  and  (2.51)  we  come  to 


the  interpretation  that  Vs  is  the  rate  of  dissipation  and  (y/4)DII/Dt 
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I'  > 0 in  a second-order  fluid  unless  3 = - y [see  equation  2.50].  But 


experimental  data  of  Markovitz  (1965)  on  polymer  solutions  do  not 


justify  3 = - y.  It  would  probably  be  a good  idea  to  use  quotation 


marks  and  say  that  Vg  is  the  "dissipation.”  However,  the  second-order 


constitutive  equation  was  derived  only  as  a slow  flow  approximation  to 


real  viscoelastic  behavior,  and  for  small  rates  of  strain  III  will  be 


small  compared  to  II.  From  (2.50)  one  can  see  that  H'g  will  be  positive 


for  such  flows. 


Another  case  in  which  'ig  will  always  be  positive  is  that  of  two- 

dimensional  flows.  If  A^,  A^  and  A^  are  the  principal  values  of  i ., 
3 3 3 

then  III  = A1  + A^  + A^.  Using  the  condition  of  incompressibility 


A,  + A„  + A = 0,  it  can  be  easily-  shown  that  III  = 3 A A A , which 
12  3 12  3 


means  III  = 0 for  all  two-dimensional  flows. 


Some  remarks  of  Coleman,  Duffin  5 Mizel  (1965)  may  be  of  interest 


here:  " ...  Truesdell  suggested  that  T. .U.  . > 0 be  a restriction  not 

ij  i,J  ~ 


on  constitutive  equation  but  rather  on  allowable  motions.  Coleman  has 


shown  that  thermodynamics  requires  only  that  T. .U.  . be  > 0 for  certain 

1 ' ij  i,J  “ 


motions  of  incompressible  fluids;  among  these  are  the  substantially 


stagnant  motions.  Coleman  has  taken  the  point  of  view  that  if  j 


be  < 0 in  one  of  these  special  motions,  then  the  constitutive  equation 


should  be  rejected.  Recently  Noll  (private  communication)  has  shown 


that,  for  a fluid  obeying  [the  second-order  constitutive  equation] 


exactly,  T. .U.  . is  > 0 for  substantially  stagnant  motions  if  and  only 

i,3  - ■ 7 


if  both  p > 0 and  3 = - y." 


We  would  like  to  make  two  comments  about  the  above  remarks  of 


Coleman,  Duffin  5 Mizel.  Firstly,  from  our  findings  here  it  seems  that 
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T. .U.  . > 0 need  not  be  a restriction  on  either  the  constitutive 

ij  i.J 

equation  or  the  al’O'-.'blt  motions,  only  the  part  of  T„IL  j which  docs 
not  go  to  elastic  energy  need  be  positive;  however  the  second-order 
fluid  fails  to  satisfy  his  unless  8 = - y.  Secondly,  we  do  not  know 
what  Coleman,  Duff’"  § Mizel  mean  by  "substantially  stagnant  motions," 
but  the  private  communication  of  Noll  referred  to  in  their  paper  agrees 
with  our  findings  if  it  means  a "motion  with  DII/Dt  = 0." 

Before  closing  the  section  it  may  be  interesting  to  try  to  separate 
the  viscous  and  elastic  parts  of  the  second-order  constitutive  equation. 
According  to  Truesdell  f)  Toupin  (1960),  p.  571,  the  stresses  in  a 
viscoelastic  fluid  can  be  divided  into  an  elastic  and  a dissipative 


part 


T . = a + T? 

ij  ij 


The  dissipative  or  non- recoverable  part  is  given  by 

V = T°.  U. 

where  T is  the  rate  of  dissipation,  and  the  elastic  or  recoverable  part 
is  found  from  the  rule  that  the  work  of  elastic  stresses  must  equal  the 
rate  of  change  of  potential  energy: 


P(PE)  _ XE  „ 

Dt  ij  i,  j 


It  is  not  possible  to  do  this  for  an  Oldroyd  or  Maxwell  fluid  in 
which,  due  to  the  presence  of  the  time  derivative  of  the  stresses, 
the  stresses  cannot  be  expressed  explicitly  in  the  form  T. . = T. . 
(kinematics) . ^ 


MB 


For  a second-order  fluid  we  have  shown  that  T .U.  . = (Y/4)DII/Dt.  and 

ij  i,j 

a little  reflection  on  its  derivation  shows  that 


_ DA.  . 

^ ■ * TE1 


'ij  * 8AikAkj  * '<(\kYj  * Ajkuk,i> 


2.12  Elastic  Energv  of  Walters  Liquid  with  Short  Memory 


The  constitutive  equation  for  the  Walters  liquid  B for  short 


memory  was 


Txj  = - Pg1:i  + 2pE1^  - 2K  — i 
6 ot 


^ + UmEij  - U1  Emj  - Uj  E11 

ot  3t  ,m  ,m 


The  deformation  work  is  therefore 


T1JU.  . = 2uE1JU.  . - 2K 


i,j  3t 


J U.  .+  UmE1J’u.  (U1  Emj’  + 


i,j  ,m  i,j  ,m 


Uj  Eim)U. 


•=  2yE1-*E.  . - 2K  % (E1^  E . .)  + h Um(E1;iE.  .)  -U1  Em;iE.  .-lP  EimE . . 

ij  3t  ijJ  A ,m  ,m  ij  ,m  j.j 


= 2pl I - 2K|i  - 2 III 


where  the  irvariants  of  the  strain  rate  tensor  E1  are  II  = E^E.  . and 

J iJ 

III  = E1Em^E...  (Note  that  in  the  previous  section  II  and  III  referred 
m ij  ^ 

to  the  invariants  of  A. . = 2E..1.  The  deformation  work  for  the  Walters 

ij  ij 

liquid  B with  short  memory  can  therefore  be  written  as 
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where 


TljU. 


D 1 1 
Dt 


V = 2 u 1 1 + 4 K III 


(2.53) 

(2.54) 


These  expressions  look  very  similar  to  the  corresponding  equations  for 

a second-order  fluid,  equations  (2.49)  and  (2.50).  In  a general 

motion  III  can  be  of  either  sign  so  that  the  dissipation  in  a Walters 

liquid  with  short  memory  cannot  be  proved  to  be  positive.  However, 

short  memory  is  the  same  as  slow  motion,  for  which  III  <<  II  and 

therefore  . > 0. 

W 


2.13  Elastic  Energy  of  Maxwell  Liquid 

Consider  first  the  Maxwell  constitutive  equation  referred  to  a 
convected  frame  of  reference 

3a . . 

* »,  if  ■ (2-55> 

where  a..(C,t)  is  the  extra-stress  and  e..(£,t)  is  the  rate  of  strain, 
xj  - ij  - 

both  referred  to  convective  coordinates;  3/3t  denotes  differentiation 
with  the  convective  coordinates  constant.  It  is  well-known  that  the 
physical  model  of  the  Maxwell  constitutive  equation  is  a spring  and  a 
dashpot  in  series  (Fig.  2.2).  Let  e denote  the  elongation  of  the  end- 
points, S is  the  force,  G is  the  spring  constant,  n is  the  dashpot 
damping  coefficient,  E = e is  the  rate  of  elongation  (velocity)  of  the 
endpoint.  Then  the  elongations  of  the  spring  and  dashpot  are  given  by 

e_  = S/G  and  e = S/n,  so  that  the  total  elongation  is 

on  " 


G 

S o—i\/\/\f 


Figure  2.2 


Model 
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e = e„  + e_  = tt  + 


/ I 


I 


assuming  that  = 0 at  t = 0.  Differentiating  the  above  equation  and 
writing  X = n/G  we  get 

s + x — = nE 

l 

which  looks  similar  to  the  Maxwell  constitutive  equation  (2.55)  if 
force  and  velocity  of  the  end-points  in  the  model  are  regarded  as 
analogous  to  the  stresses  and  strain  rates  in  a viscoelastic  material. 
Let  us  find  the  rate  of  dissipation  in  the  dashpot  and  the  rate  of 
increase  of  elastic  potential  energy  PE  in  the  spring: 

Rate  of  Dissipation  = (Force  on  Piston)  (Velocity  of  Piston) 

2 

c * S 

= S e = — 

n n 


3PE 

at 


= (Force  on  Spring)  (Rate  of  Elongation  of  Spring) 


= S e = S(e-e  ) = S(E  - £) 

u n n 


3PE  s s 

Rate  of  Total  Work  = Rate  of  Dissipation  + = — + S(E  - — ) 

d l n n 


SE 


With  the  above  analogy  in  mind,  we  suggest  that  the  rate  of  dissipation 
and  the  rate  of  increase  of  PE  for  the  Maxwell  fluid  (2.55)  are 


Rate  of  Dissipation  = 


a . .0 
il 

2u 


a.  . 


9£E_(§,t)  ij  { _ _LL) 

at  1 ij  2y  > 


(2.56) 

(2.57) 
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Leaving  aside  the  analogy,  the  above  suggestion  must  be  correct  because 
by  taking  a process  sinusoidal  in  time  given  by  a sin  ut,  it  can  be 
shown  that  the  time  integral  / a^o1-1  dt/2y  always  increases  with  time, 

whereas  PE  = a1-*  (e„  - (ck  j/2jj))dt  may  be  positive  or  negative, 

the  integral  over  a complete  period  being  zero.  This  is  shown 
rigorously  in  the  next  section  for  a fluid  more  general  than  the  Maxwell 
fluid. 

Now  referring  (2.55),  (2.56)  and  (2.57)  to  fixed  coordinates  we 
conclude  that  for  the  two  kinds  of  Maxwell  fluids 


<SS . . 

S.  . + A -ril  = 2pE. . 

lj  i 6t  ij 


(2.58) 


S1*  ♦ A 

i 


= 2pE 


(2.59) 


the  dissipation  and  PE  are  given  by 


S.  .S1J 

ii 


(2.60) 


P(PE)  _ sijrE 

Dt  1 ij 


(2.61) 


In  (2.58)  and  (2.59),  5/5t  denotes  the  convective  derivative  of  Oldroyd, 
and  it  has  been  replaced  by  D/Dt  = 3/3t  + U^(  ) ^ in  (2.61)  because 
when  operating  on  a scalar  <5/6 1 = D/Dt.  The  deformation  work  is 


sij.E  - ? + £IpE) 

b ij  ¥M  Dt 


(2.62) 


as  is  easily  seen  by  adding  (2.60)  and  (2.61).  Equations  (2.58) 
through  (2.62)  are  the  results  of  this  section. 
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2.14  Elastic  Energy  of  Qldrovd  Liauid 


Consider  the  Oldroyd  equation  referred  to  convected  coordinates. 


3o.  . 

o.  . + X 
1J  i at 


3e.  . 

= 2u(e  + X — 

lj  2 


(2.63) 


The  physical  model  of  this  constitutive  equation  [see  Burgers  (1935)] 
is  shown  in  Fig.  2.3.  Let  suffixes  1,  2 and  3 respectively  refer  to 
the  series  dashpot,  the  parallel  dashpot  and  the  spring,  as  shown. 
Then  the  force  S at  the  end  is  given  by 


S = S + S = n e + G e 
2 3 2 2 3 3 


Using  the  conditions  e=e  +e,e  = e and  S = n e , the  above 

12-23  11 

becomes 

• n2  G3  « 

S = n (e-e  ) + G (e-e  ) = n e - — S + G e - — / Sdt 


2 n. 


3 n o' 


which  can  be  rewritten  as 


t 

/ Sdt 
% ' 


Vn2 


n2  • 

s = n (e  + rr-  e) 
1 G3 


or,  differentiating,  as 


s + y i£ 
s lit 


rc  i 3E. 

n (E  + X rr) 
1 2 ot 


(2.64) 


where  we  have  put  E = e and 


and 


Vn2 


(2.65) 


(2.66) 


The  model  equation  (2.64)  is  analogous  to  the  Oldroyd  constitutive 
equation  (2.63)  in  convected  coordinates. 


Now  consider  the  rate  of  dissipation  in  the  two  dashpots  of  the 


model : 


Rate  of  dissipation  = Se  +Se  =Se  + S (e-e  ) 
r 12  2 12  1 


= (S-S  ) e , + S e = (S-S  ) - + S e 

2 1 2 v 2 n 2 

2 

= — + S (e  - — ) 
n 2 n 

i i 


Since  S = n e = n (e-e  ) =n  (e  - S/n  ),  the  above  becomes 
2 2 2 2 1 2 1 


Rate  of  Dissipation 


r—  + n (E  - r—) 

"i  2 \ 


Now  consider  the  rate  of  increase  of  PE  in  the  spring: 


= S e = (S-S  ) (e-e  ) = [S  - n (e  - ~ ] (e 
33  21  2 H , 


c c 2 

. S(E  - 4-)  - n (E  - -f) 

nl  2 "l 


The  constant  n may  be  expressed  in  terms  of  X and  X by  dividing 
2 12 


(2.46)  and  (2.47),  which  yields 


X -X 
1 2 


(2.67) 


The  expressions  for  dissipation  and  elastic  energy  for  the  model  then 


become 


We  are  not  saying  anything  about  the  significance  of  A and  B now.  To 
find  that  out,  assume  that  a material  particle  undergoes  a sinusoidal 
process  in  its  path,  represented  by 
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e . . (£ ,t)  = F.  . (£)  cos  wt 

lj  - ' ij  *• ~ 

where  F(£)  remains  constant  along  the  path  of  the  particle.  Now  the 
part  A is  non-negative,  and  hence  its  time  integral  will  always  grow 
with  time.  If  we  can  show  that  the  time  integral  of  B over  a complete 
period  is  zero,  then  we  can  conclude  with  some  assurance  that  A is  the 
rate  of  dissipation  of  energy  and  B is  the  rate  of  change  of  potential 
energy. 

Substitution  of  the  process  = F^  cosoit  into  the  constitutive 
equation  (2.63)  results  in  a linear  nonhomogeneous  equation  for  , 
whose  particular  solution  is 


a.  . = 

ij 


2F..uu>(A  -A  ) 2F.  . y (1  + A.A  “ 

— zl — sinwt  + —j cosuit 

1 + A^w  1 ♦ X u> 


and  whose  homogeneous  solution  is  cr^  = (const.)  exp(-t/A^);  this  of 
course  dies  out  very  quickly  (it  will  be  zero  exactly  if  the  motion 


started  from  rest)  and  therefore  will  not  be  considered.  With  the 
stresses  and  strain  rates  thus  determined,  a straightforward  integration 
shows  that 


Bdt  = 0 


We  therefore  conclude  that  A is  the  rate  of  dissipation  and  B is  the 
rate  of  change  of  PE  for  the  Oldroyd  fluid  in  convected  coordinates. 

Now  referring  (2.65),  (2.70)  and  (2.71)  to  fixed  coordinates  we 
conclude  that  for  the  Oldroyd  fluids 


6 S.  . 

S.  . + A -rpJ- 

ij  i 6t 


6E.  . 


(2.72) 
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(2.73) 

(2.74) 

(2.75) 


The  deformation  work  is 


S.  .E1-* 
ij 


D(PE) 

Dt 


(2.76) 


as  is  easily  verified  by  adding  (2.74)  and  (2.75).  Equations  (2.72) 
through  (2.76)  are  the  results  of  this  section. 


CHAPTER  3 

STABILITY  OF  PLANE  COLETTE  FLOW 


In  this  chapter  the  stability  of  a plane  Couette  flow  of  the 
second-order  fluid  will  be  studied  by  the  variational  technique  of  the 
energy  method.  It  is  first  necessary  to  derive  the  energy  equation. 


3.1  Equation  of  Global  Disturbance  Energy 

The  starting  point  of  the  energy  method  is  the  equation  of  dis- 
turbance mechar  ical  energy,  integrated  throughout  the  flow  field.  To 
derive  this,  let  the  velocity,  pressure  and  stress  tensor  of  a disturbed 
motion  be  broken  up  into  a basic  steady  motion  and  a disturbance  (not 


necessarily  small): 


V.  = U.  u. 
111 


IT  = P + p 


(3.1) 


t.  . = T.  . + t.  . 

The  disturbances  vanish  on  the  boundaries  and  are  either  periodic  or 
invariant  in  the  directions  in  which  the  flow  is  unbounded.  For 
example,  for  a plane  Couette  or  Poiseuille  flow,  the  disturbances  are 
either  periodic  or  invariant  in  the  direction  parallel  to  the  walls. 

Since  both  the  total  and  the  basic  motion  must  satisfy  the 
equations  of  motion,  we  have 


O V . ■. 

— -1-  + V.V.  . = - T..  . 

at  j i,j  p ij.j 


u.u.  . = - t.  . 

j i,j  p ij.i 


(3.2) 
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Subtracting  the  above  equations  tve  get  the  equation  of  motion  of  the 
disturbance 


3u. 

i 

3t 


— t..  . - u.U.  . - U.u.  . - u.u. 
P 13,3  J i,3  3 i,3  3 


Multiplying  the  above  by  u^,  summing,  and  integrating  over  a volume  V 
fixed  in  space,  we  get  the  energy  equation 


/“  X 


3u. 

l 

at 


dV 


— fu.t. . . dV  - fu.u.U.  . dV  - fu.U.u.  . dV 
P J i iJ,J  J i 3 i,J  1 i 3 1,3 


J u . u . u . . dV 
; i 3 1,3 


Let  us  choose  the  volume  V so  that  it  coincides  with  the  walls  and 
extends  to  a length  equal  to  an  integral  multiple  of  the  wavelength  in 
the  direction  of  periodicity  (Fig.-  3.1).  Then  with  the  help  of  the 
equation  of  continuity  (IL  ^ = 0;  u^  ^ = 0) , and  Gauss's  theorem 
relating  surface  and  volume  integrals,  the  terms  on  the  right  side  of 
the  energy  equation  can  be  reduced  as  follows: 


fu.t. . .dV  = /( u.t. .)  .dV  - ft. .u.  .dV  = |(t. . u.)dS.  - ft. -u.  . dV 

> l 13 ,3  1 1 13  ,3  1 13  1,3  J 13  1 3 J 13  1,3 


fu.U.u.  -dV  = k f U.(u2)  .dV  = % /( u2U.)  .dV  = h { u2U.dS. 

1 1 3 i*3  J 3 1 ,3  1 i 3 ,3  y 1 3 3 

2 2 2 
fu.u.u.  - dV  = / u . (u. ) . dV  = '2  /(u.u.)  . dV  = ^ / u.  u . dS  . 

; 1 3 1,3  J 3 1 ,3  J 1 3 ,3  r 1 3 3 


The  transport  terms  (surface  integrals)  in  the  above  are  zero  because 
the  integrands  are  zero  on  those  portions  of  V coinciding  with  the  walls, 
and  their  contributions  cancel  each  other  on  the  open  portions  of  V 
(see  Fig.  3.1  again).  The  equation  of  energy  therefore  reduces  to 


r 
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2 

% f u . d V = - — / t..u.  .dV-  / u.u.U.  •( 
dt  1 l P J ij  i,J  1 l J 1,3 


[At  this  point  note  that  the  Newtonian  equation  (1.1)  is  obtained  if 

t^j  = - P^j  + V (u^  + Uj  ^)  is  substituted  in  the  above  equation.] 

Introducing  the  strain  rate  tensor  of  the  basic  flow  A. . = U.  . ■*■  u. 

ij  i.J 

the  above  becomes 


7-4-/  u.dV  = - % f u.u.A.  .dV  - / t..u.  .dV 

2 dt  ^ l 2 ; l j ij  J ij  i,j 


(3.3) 


This  is  the  global  energy  equation  for  the  disturbance,  and  is  attributed 
to  Reynolds  (1895)  and  Orr  (1907).  It  signifies  that  the  rate  of 
change  of  kinetic  energy  of  the  disturbance  in  the  entire  field  equals 
the  energy  extracted  from  the  basic  flow  by  means  of  the  Reynolds 
stresses  - minus  the  "dissipation"  of  the  disturbance  kinetic 

energy  into  heat.  We  have  jut  the  word  "dissipation"  within  quotation 
marks  here  because  from  the  discussions  of  Sections  2.10-2.14  it  will 
be  apparent  that  for  viscoelastic  fluids  which  can  store  energy,  this 
term  represents  the  dissipation  rate  as  well  as  the  rate  of  change  of 
elastic  potential  energy  of  the  disturbance. 

As  an  aside  it  may  be  pointed  out  that  for  the  application  of  the 
energy  method  to  Newtonian  fluids  one  poses  the  following  problem: 
consider  the  class  of  disturbances  having  instantaneously  stationary 
global  kinetic  energy;  what  is  the  equation  obeyed  by  that  member  which 
can  exist  with  the  largest  viscosity  (other  things  being  equal)?  Thus, 
one  puts  the  left  hand  side  of  (3.3)  to  zero,  introduces  t^  = - p6 + 

p(u.  . ♦ u.  .)  and  maximizes  u by  the  calculus  of  variations.  The 

* > J J I1 


83 


Euler  equation  corresponding  to  this  variational  problem  represents  the 
most  efficient  or  most  dangerous  disturbance,  one  which  can  extract 
energy  from  the  basic  flow  most  efficiently,  losing  as  little  as 
possible  to  the  viscous  dissipation. 

3.2  Equation  of  Global  Disturbance  Energy  for  Second-Order  Fluid 

It  was  found  that  due  to  the  presence  of  the  time  derivative  in  the 
expression  for  t^  for  the  second-order  fluid  [see  equation  (3.5)  later], 

the  energy  method  as  formulated  for  a purely  viscous  fluid  is  not  appli- 

* 

cable  to  a second-order  fluid.  As  a way  around  this  difficulty,  Lumley 
in  a private  communication  suggested  that  the  time  derivative  can  be 
avoided  if,  instead  of  considering  disturbances  having  instantaneously 
stationary  kinetic  energy,  one  considers  those  having  instantaneously 
stationary  tc.  al  (kinetic  + elastic)  energy. 

Moreover,  the  consideration  of  the  total  energy,  instead  of  simply 
the  kinetic  energy,  seems  to  make  better  sense  in  a viscoelastic  fluid. 
The  disturbance  in  a viscoelastic  liquid  is  at  all  times  interchanging 
energy  between  its  kinetic  and  elastic  (potential)  modes.  At  one  moment 
it  may  h.ave  very  little  kinetic  and  a good  deal  of  potential  energy, 
while  at  another  moment  it  may  have  the  opposite.  Therefore  it  is  the 
total  (kinetic  + potential)  energy  that  must  be  significant  for  a dis- 
turbance in  a viscoelastic  fluid. 


Note  that  even  in  the  stationary  neutral  state  of  no  rate  of  change  of 
total  disturbance  kinetic  energy  in  the  flow  field,  the  3/9t  of  a 
quantity  is  not  zero.  This  is  because,  although  the  integral  over  the 
entire  field  is  zero,  the  rate  of  change  of  kinetic  energy  at  every 
point  is  not;  at  some  points  it  is  increasing,  while  at  others  it  is 
decreasing.  This  must  happen  due  to  the  nonlinearity  of  the  dis- 
turbances in  the  energy  method  (Lumley,  1971b). 
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Recall  from  Section  2.11  that  the  elastic  potential  energy  of  a 
disturbance  in  a second-order  fluid  is  yIi'/4  where  II ' is  the  second 

I 

invariant  for  the  strain  rate  tensor  of  the  disturbance  motion:  II  = 

a. .a. . where  a. . = u.  . + u.  ..  With  this  in  mind  let  us  rewrite  (3.3) 

ij  ij  i,J  l»i 


as 


2 

u. 


& / T dv  * } / 5T  dv  ' - I /“iV«dV  - i 5T  > 


Y DII 


dV 


Inserting  DII  /Dt  = 3 1 1 /3t  + V^II  ^ and  using  continuity  and  Gauss's 
theorem,  the  second  integral  on  the  left  hand  side  can  be  rewritten  as 
follows: 

t i 

/ BT  dv  * / W-  dv  " /(VI'>,kdv  - ^ / n'dv  . fvkn'dsk 


The  surface  integral  in  the  above  equation  vanishes  because  of  boundary 

conditions,  so  that  the  energy  equation  becomes 
2 

Y Dll' 


j pu.  i 

4?  / Hr  + 4 II  )dV  = - / u.u.A.  .dV  - /( t.  .u.  . - 5-  ~ )dV 

dt  J K 2 4 1 2 J l j ij  J 1 ij  i,j  4 Dt 


(3.4) 


That  is,  the  rate  of  change  of  the  sum  of  the  disturbance  kinetic  and 
elastic  energies  equals  the  total  production  (through  Reynolds  stresses) 
minus  the  total  dissipation. 

An  expression  for  the  stress  change  t^  which  occurs  in  the 

equation  of  energy  will  now  be  found  for  a second-order  fluid.  Let 
(2) 

and  A. . be  respectively  the  first  and  second  Rivlin-Ericksen  tensors 
1J  (2) 

of  the  basic  motion,  and  let  and  be  those  for  the  total  motion. 

Thus , 


A. . = U.  . + U.  . 

ij  ill  J.i 


where  a. . = u.  . + u.  ..  Substitution  of  the  above  expression  for 

ij  J.i  * 

stress  change  in  (3.4)  results  in  the  equation  of  global  disturbance 
mechanical  energy  for  the  second-order  fluid. 


I 

3.3  The  Extremization  Problem 

Consider  the  following  extremization  problem:  out  of  the  class  of 

disturbances  having  instantaneously  stationary  total  (kinetic  + elastic) 
energy,  what  is  the  equation  obeyed  by  that  member  which  can  exist  with 
the  largest  viscosity  p (other  things  being  equal),  i.e.,  is  most 
dangerous  or  most  capable  of  causing  instability? 


* 
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Substituting  the  stress  change  (3.5)  into  the  energy  equation 
(3.4),  we  get  for  the  instantaneous  stationarity  of  disturbance  energy 

- | /WijdV  = “ NUi.jdV  * 8 haikAkj  * Aikakj  * aikakj>ui.jdV 

|ukAij,k  * aikUk,j  * Aik\j  * aikuk,j  (3'6) 


* ajk^k,i  + ^k^.i  + ajkuk,i 


u.  .dV 
i.J 


The  viscosity  y will  now  be  extremized  in  the  above  equation.  The 
question  is,  what  happens  to  the  coefficients  0 and  y as  p is  varied? 
If  we  imagine  that  the  viscosity  of  a liquid  containing  polymer  mole- 
cules is  gradually  raised,  we  expect  that  the  polymer  molecules  will 
become  more  and  more  effective  and  therefore  will  increase  the  non- 
Newtonian  parameters  0 and  y.  It  is  simplest  to  assume  that  the 
variations  are  proportional: 


60  = Kg6y 

6y  = K 6u 
Y 


(3.7) 


where  Ka  and  K are  constants.  It  will  be  shown  in  the  next  section 
6 Y 

that  the  results  of  the  present  work  can  be  justified  through  another 
type  of  extremization. 

The  viscosity  y in  (3.6)  is  extremized  by  the  method  of  calculus 
of  variations;  that  is,  take  the  variation  of  each  term  of  the  equation 
and  set  6y  = 0 for  the  extremum.  During  the  computation  of  the 
variation  of  the  different  terms  Gauss's  theorem,  boundary  conditions 
and  the  equation  of  continuity  are  used  several  times.  The  algebra  is 
long  but  perfectly  straightforward.  The  result  is 


87 


2yaij  + ^2aikAkj  + 2Aikakj  + 3aikakj}  + Y(ukAij,k  * aikAkj 

K “k./ik  + \,i\j  + ^ikVj.j  + YUk,jAkj,i  6uidV  = 0 

J 


Following  the  procedure  of  Lumley  (1971b),  note  that  in  the  above 
equation  the  variations  are  over  incompressible  motions  so  that  the 
6u^  are  not  independent  but  are  related  by 


6u.  . = 0 


(3.9) 


We  may  specify  6u^  and  6u^  independently,  for  example,  and  6u^  is  then 
given.  Denoting  the  quantity  within  the  braces  { } in  (3.8)  by  X^,  we 


can  write 


fx.6u.dV  * f(X  Su  + X 6u  + X 6u  )dV  = 0 
J 1 1 J 1 1 2 2 3 3 


If  we  set  X = - $ , where  <ji  is  any  arbitrary  function,  then  integra- 

3 , 3 

tion  by  parts  (using  boundary  conditions)  and  use  of  (3.9)  gives 

0 = fX.6u.dV  = f(X  6u  + X 6u  + <j>6u  )dV 

‘ 1 1 11  2 2 3,3 

* / (X  6 u +X  6u  -<t>[6u  +6u  ])dV 

1 1 2 2 1,1  2,2 

= /(X  6u  +X  6u  +$  6u  -*-4>  6u  )dV 

112  2,11,22 

= /((X  +$  )6u  + (X  +♦  )6u  ]dV 

' 1 ,1  1 2 ,2  2 

Since  6u  and  6u  are  independent  we  must  have  X.  = - $ ..  From  the 

12  ‘ 1 ,1 

definition  of  X^  we  therefore  get 


f 
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pAijUj  * -v 


2paij+  e(2aikAkj  + 2Aikakj  + 3aikakj5  + Y(ukAij,k 


♦ a..  A,  . + 2A..a,  . + u,  .A.,  + u.  . A.  . + 3a.,  a,  . ) 
lk  kj  lk  kj  k,j  lk  k,i  kj  . lk  kj' 


,j  “ YUk» jAkj 


(3.10) 

This  is  the  Euler- Lagrange  equation  corresponding  to  our  variational 
problem,  and  is  satisfied  by  the  most  dangerous  disturbance  in  a second- 
order  fluid  obeying  assumptions  (3.7).  The  same  equation  is  also 
obtained  if  one  uses  the  technique  of  Lagrange  multiplier,  i^i  which 
case  p would  be  identified  as  the  Lagrange  multiplier  satisfying  the 
constraint  of  incompressibility.  Also  note  that  except  for  the  last 
term  which  is  not  a derivative  with  respect  to  x^  , equation  (3.10)  looks 
very  similar  to  an  equation  of  motion,  in  which  case  <j>  appears  like  a 
"pressure."  Lumley  (1971b)  for  the  Newtonian  case  uses  the  expression 
"$  is  the  'pressure'  which  must  be  applied  to  the  fluid  in  order  to 
satisfy  incompressibility." 


3.4  Extremizing  Concentration  Instead  of  Viscosity 

In  the  last  section  u was  extremized  and  the  variations  of  6 and  y 
were  taken  to  be  proportional  to  that  of  y.  It  would  be  interesting  to 
see  what  would  happen  if  we  took  a linear  dependence  of  y,  B and  y on 
the  concentration  C of  polymer 

V = V * wC 

0 

8 = Tc 

Y = yC  (3.7)' 


89 


where  y,8,y  are  constants  of  proportionality,  and  then  maximized  the 
concentration.  Physically  this  variational  problem  would  mean  that  as 
more  and  more  polymer  molecules  are  being  added  to  a turbulent  Couette 
flow  y,  8 and  y are  all  increasing  in  magnitude;  a concentration  will 
eventually  be  reached  at  which  the  disturbances  just  begin  to  vanish. 

With  the  assumed  dependence  of  y,  8 and  y on  C,  if  the  stress  change 
of  (3.5)  is  substituted  in  the  energy  equation  (3.4),  we  get  for  the 
disturbances  having  instantaneous  stationarity  of  energy 


- f/“iuj''ij<lv*,‘0/aijui,jdv*c/[lraij  * staikAjk  * Aikajk*aikajk)*i'(ukAij,k 


* aikukj  * Aik\j  * aik\,j  * ajkuk,i  * Ajkuk,i*ajkuk,i>'uijdV 


Taking  the  variation  of  each  term  of  the  above  equation  and  setting 
6C  = 0 for  maximum  C,  we  have  found  that  equation  (3.10)  of  the  last 

i 

section  results.  This  is  not  surprising  since  both  (3.7)  and  (3.7) 
amount  to  the  fact  that  the  variations  of  8 and  y are  proportional  to 
that  of  y. 

Note  also  that  if  we  had  regarded  8 and  y as  constants  as  y was 
varied,  we  would  again  have  obtained  the  same  result,  that  is  (3.10). 
This  is  because  at  the  critical  state  we  set  5y  = 0,  which  makes  66  = 

6y  = 0 also. 

Thus,  the  results  of  this  chapter  can  be  justified  from  three 
points  of  view. 

3.5  Rectilinear  Couette  Flow 

Let  us  consider  the  special  case  when  the  basic  flow  is  in  the  form 
of  a plane  Couette  flow  of  shear  rate  P (see  Fig.  3.2): 
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u = rx 

1 z 

u = u = o 

2 3 


(3.11) 


For  this  basic  motion,  it  follows  from  the  constitutive  equation  of  the 
second-order  fluid  that  the  stress  tensor  is 


[T]  = 


-P+(g+2y)r 


An  examination  of  the  equation  of  motion  then  shows  that 


P = constant  throughout  the  flow 


For  the  basic  motion  of  (3.11),  the  three  equations  in  (3.10)  simplify 


pTu  = -<i>  + 2\i  u + B[2r(3u  + u + u ) + 3(a  a,  .)  . ] 

2 1, j j 1,12  2,11  2 , j j Ik 


+ v[r(3u  + 4u  + 2u  ) + 3(a  a.  .)  .] 
2, j j 1,12  2,11  lk  kJ  » J 


pTu  = -<j>  + 2p  u + S[2r(3u  +u  +u  ) + 3 (a  a,  .).] 

1 ,2  2 , j j 2,12  1,22  l,jj  2k  kJ  ,J 


♦ Y[r(3u  + 4u  + 2u  ) + 3(a  a,  .)  .] 
1 i,jj  2,12  i,22y  2k  kr  ,y 


0 = - + 2u  u + B[2T(2u  + u + u ) + 3(a  a,  .)  .] 

,3  3, j j 3,12  1,23  2,13  3k  kJ 


+ y[2r (u  + u + u ) + 3(a  a.  ) .]  (3.12) 

3,12  1,23  2,13  3k  kJ 


Two  special  cases  will  now  be  considered,  namely  that  of  the  two- 
dimensional  disturbances  and  that  of  the  "longitudinal  rolls." 
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3.6  Formulation  of  the  Two-Dimensional  Disturbance  Problem 

Consider  only  two-dimensional  disturbances  for  which  u^  = 0 and 

( ) = 0,  so  that  the  third  equation  of  [3.12)  is  absent.  A stream 

» 3 

function  can  be  defined  through 

u = 4 

1 .2 


U = -U) 

2 .1 


Let  us  first  evaluate  the  terms  (a.,  a,  .)  . which  occur  in  [3.12) 

IK  KJ  , J 


[aikVj> 


u +u 

1,2  2,1 


u +u 

2,1  1,2 


u +u 

1,2  2,1 


u +u 

2,1  1,2 


4»  -4 

,22  ,11 


4>  -4, 

,22  ,11 


4,  -4, 

,22  ,11 


4/  -4, 

,22  ,n 


441  4)  + (4  -41  ) 

,12  ,12  ,22  ,11 


44)  4 1 +(4>  -<p  ) 

,12  ,12  ,22  ,11 


For  i = 1 and  2 this  gives 

(a.,  ^ J 4 = (a..  a,  .) 


a ) = 

[44/  4 + 

(4, 

' * , 

lk  kl  ,1 

,12  ,12 

,22 

,11 

a ) 

[44)  4)  + 

(4, 

- 4* 

2k  k2  ,2 

C\J 

A 

CM 

,22 

In  terms  of  the  stream  function  the  first  two  equations  of  (3.12) 
therefore  become 
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-pTijj  = -4>  + 2p4)  +0{2T  C3^ 

,1  ,!  ,2jj 

+ y(T(-34)  +44)  -2<p 

,ljj  ,122 

pT4)  = -<()  -2yip  +0{2r(-34) 

,2  ,2  ,ljj 

+ y{r(34)  - 4 4)  +24) 

> 2 j j ,H2 


-4)  -4)  )+3[44)  4»  +('|)  -4>  ) ] 

122  ,111  ,ljj  ,12  ,12  ,22  ,11  ,1 


) + 3 [44)  4)  + (4>  -4)  ) ] } 

111  ,12  ,12  ,22  ,11  ,1 


+4)  +4)  ) + 3[44)  4)  +(4)  -4)  ) ] 

112  ,222  ,2jj  ,12  ,12  ,22  ,11  ,2 


) + 3 [44)  tp  +(4)  ~4»  3 ] } (3.14) 

222  ,12  ,12  ,22  ,11  ,2 


Elimination  of  the  "pressure"  $ by  cross-differentiation  and  subtraction 
between  the  above  equations  results  in  the  cancellation  of  all  the  non- 
linear terms  and  the  terms  containing  3,  giving 

pT4)  + u4)  ••••  - 2yF\p  = 0 (3.15) 

,12  ,HJ1  , 12j  j 

The  final  equation  is  therefore  linear,  although  the  problem  is  not. 

The  boundary  conditions  are  u =u  = 0 at  x = ± 6 (Fig.  3.2).  In 

1 2 2 

terms  of  the  stream  function  the  boundary  conditions  read  4)  =4)  =0 

,1  ,2 

at  x = ± <5.  But  4)  = 0 at  x = ± 6 implies  4*  = constant  at  x = ± 5: 

2 ,12  2 

for  convenience  we  take  the  constant  to  be  zero.  The  boundary  con- 
ditions ‘are  therefore 

4)=4>  =0  at  x=±6  (3.16) 

,2  2 

Equation  (3.15)  with  boundary  conditions  (3.16)  define  a characteristic 
value  problem  which  will  now  be  solved  following  the  procedure  of  Orr 
(1907)  who  had  equation  (3.15)  with  y = 0. 
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3.7  Solution  of  the  Two-Dimensional  Disturbance  Problem 

To  solve  the  characteristic  value  problem  (3. 15) - (3. 16) , note  that 
since  the  disturbances  are  periodic  in  the  x^-direction 


. » , . tax, 

<l>  = ^ (x  ) e 1 
2 


(3.17) 


where  the  wavenumber  a is  real.  Substituting  the  above  in  (3.15)  and 
expressing  in  terms  of  non-dimensional  distance  and  wavenumber 


; = 


k = a5 


(3.18) 


we  get  (writing  D = d/d;) 


2 2 2-  - 2 2. 

(D  -k  ) tj>  + ikRD^  - 2ikGD(D  -k  )<j;  = 0 


(3.19) 


where  we  have  also  introduced  the  non-dimensional  numbers 


R = Reynolds  Number  = 


2 

p6  r 


2 2 2 

P U 6 2 

I = (5  ) 


(3.20) 


r ...  . . ..  n . yr  time  scale  of  fluid 

G = elasticity  Parameter  = J : -r— r, 

J y time  scale  of  flow 


In  the  above  definition  of  Reynolds  number,  is  called  the  "shear 
velocity"  and  is  defined  as 


2 

u 

T 


wall  shear  stress  _ jiT 
density  p 


and  S+  = pSu^/y  is  the  non-dimensional  half-width  of  the  channel. 

Since  the  coefficients  of  (3.19)  are  independent  of  the  solution 
of  the  equation  can  be  expressed  as  a superposition  of  solutions  of  the 
form  4 1 = elmJ\  With  this  substitution  (3.19)  becomes 
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4 3 2 2 2 

m - 2Gkm  + 2k.  m - k(2Gk  + R)m  = 0 (3.21) 


Denoting  the  roots  of  this  equation  by  n^,  m , m3  and  m 4,  the  general 
solution  is 

im^C  im2?  im3£  im4; 

= Ae  + Be  + Ce  + De 


Because  of  the  boundary  conditions  <>  = D^  = 0atC  = ±l,  the  condition 
that  not  all  of  A,  B,  C and  D are  zero  is  the  vanishing  of  the  deter- 
minant 


ltn 


1 


-lm. 


-im„ 


lm. 


-lm. 


lm, 


-lm, 


lnij  im2  im3  im4 

me  me  me  me 

12  3 4 

-im1  -im^  -in3  -im^ 

me  me  me  me 

12  3 4 


= 0 


(3.22) 


By  expanding  the  determinant  the  above  condition  can  be  shown  to  reduce 

to 


(m  m +m  m )sin(m  -m  )sin(m  -m  ) + (m  m +m  m )sin(m  -m  )sin(m  -m  ) 
3412  12  34  2413  24  31 

+ (mm  +mm)sin(m-m)sin(m  -m)  = 0 

2 3 1 4 2 3 1 4 


(3.23) 


The  procedure  of  solving  the  characteristic  value  problem  is  now  clear: 
for  an  assumed  value  of  k and  G,  find  R so  that  the  four  roots  of  (3.21) 
satisfy  (3.23);  this  gives  a point  on  the  stability  diagram.  The  pro- 
cedure was  repeated  for  several  values  of  k and  G on  the  computer  and 
the  results  are  shown  in  Fig.  3.3.  G,  the  ratio  of  the  time  scales  of 
the  fluid  and  flow,  was  given  a value  of  G = -0.5  and  -0.05.  Since  y 
(and  hence  G)  is  negative  on  thermodynamic  grounds,  the  curve  for  G = 
0.5  shown  dotted  in  Fig.  3.3  has  no  physical  meaning. 


REYNOLDS  NUMBER 


Figure  3.3  Neutral  stability  curves  for  two-dimensional  disturbances 
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- From  the  stability  diagram  it  is  seen  that  the  critical  Reynolds 
number  corresponding  to  G = - 0.5  is  Rcr  = 57.8,  and  that  corresponding 
to  G = 0 is  Rcr  = 44.3,  which  agrees  with  Orr's  solution  of  the 
Newtonian  case.  This  shows  that  the  presence  of  elasticity  stabilizes 
the  flow. 

From  the  point  of  view  of  the  viscous  sublayer  thickness,  the 
results  are  very  interesting.  It  will  be  recalled  from  Section  1.6 
that  half  the  channel  width  at  the  critical  state  is  a good  estimate  of 
the  viscous  sublayer  thickness  of  a turbulent  flow.  Since  the  non- 
dimensional  channel  half  width  5*  equals  the  square  root  of  the  Reynolds 
number  [see  equation  (3.20)],  we  have  the  results 

6+  = /44 . 3 = 6.65  for  G = 0 (Newtonian) 

<S+  = /577a  = 7.6  for  G = - 0.5 

This  represents  an  increase  of  14.3%  in  the  sublayer  thickness  and 
agrees  with  the  experimentally  observed  fact  that  the  sublayer  thickness 
increases  during  the  Toms  phenomenon. 

3. 8 Formulation  of  Three-Dimensional  Disturbance  Problem 

Now  consider  disturbances  which  have  velocity  components  in  all 
three  directions,  but  which  are  invariant  in  the  streamwise  direction 
x^.  Experiments  of  Bakewell  5 Lumley  (1967)  indicate  that  this  is  the 
type  of  disturbance  prevalent  in  the  wall  region  of  a turbulent  flow. 

The  motion  is  of  the  nature  of  streamwise  vortex  pairs  (Fig.  3.4),  the 
planes  of  circulation  of  which  are  tipped  so  that  their  normals  lie  more 
or  less  in  the  direction  of  the  maximum  positive  strain  rate.  Due  to 
the  mean  vorticity,  these  normal,  weep  from  an  orientation  more  nearly 


*2 


(o)  COMPOSITE  REPRESENTATION  (b)  STREAMLINES  IN  A PLANE 

NORMAL  TO  MEAN  FLOW 

Figure  3.4  Eddy  pair  in  the  wall  region  of  a turbulent  flow  [after 
Bakewell  6 Lumley  (1967)]. 


w 
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perpendicular  to  the  wall,  through  the  direction  of  maximum  positive 

strain  rate,  and  finally  to  an  orientation  nearly  in  the  streamwise 

direction.  As  they  do  this,  their  efficiency  of  energy  extraction  from 

the  mean  motion  first  increases,  and  then  decreases,  resulting  in  first 

a rapid  growth  followed  by  relatively  slow  decay. 

With  ( ) =0,  equations  (3.12)  become 

» 1 

pru  = 2pu  ..  + 8[2ru  _ ♦ 3(a  a.  .)  .]  + y[3ru  ..  + 3(a  a.  0 ■) 

2 1,JJ  2,jj  lk *3  >J  2, j j lk 

pTu  = - $ + 2pu  + B[2r(u  + u ) + 3(a  a.  .)  .] 

1 ,2  2,jj  1,22  l,jjJ  "2k^/,JJ 

+ Y[r(3u  . . + 2u  ) + 3(a  a..)  - u a..  J (3.24) 

1,JJ  1,22  2k  KJ  *2 


0 = - <J>  + 2yu  + e[2ru  + 3(a  a,  .)  .] 

,3  3, j j 1,23  3k 


+ Y[2ru  + 3(a  a.  .)  . - u .a  ] 
1,23  3k  ^ 


In  the  above  equations  <p  ^ has  been  set  to  zero  because  we  are  con- 
sidering disturbances  which  cause  no  streamwise  "pressure"  gradient. 

After  elimination  of  $ by  cross-differentiation  between  the  last 
two  equations  of  (3.24),  it  was  found  that  the  nonlinear  terms  do  not 
cancel,  and  the  solution  of  the  problem  seems  formidable.  For 
simplicity  we  have  solved  the  case  when  the  disturbances  can  be  assumed 
small,  so  that  the  nonlinear  terms  in  the  disturbance  quantities  can  be 
neglected.  Equations  (3.24)  then  simplify  to 

2 2 

pru  = 2p7  u + r f 23  + 3y)7  u (3.25a. 

2 1 2 


A 
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2 2 

pTu  = + 2yV  u + T(20  + 3y)V  u + 2r(S+y)u  (3.25b) 

1 .2  2 1 1.22 


0 = -<t>  + 2yV  u + 2r(g+y)u 

,3  3 1,23 


where  V * ( ) + ( ) . The  set  of  equations  (3.25)  is  to  be 

,22  ,33 

solved  subject  tou  =u  =u  = 0 at  the  boundaries. 

12  3 

Since  the  solutions  are  periodic  in  the  x^-direction,  introduce 

u = u(x  ) cosax 
1 2 3 


u * v(x  ) cosax, 


u = w(x  )sinax 
3 2 3 


(3.26) 


4>  = 6(x  ) cosax 

2 3 


where  a is  the  wavenumber.  The  sine  with  u and  cosine  with  u are 

3 2 

necessary  in  order  to  be  consistent  with  the  equation  of  continuity 


^ + u^  =0.  Equations  (3.25)  and  the  continuity  then  become 


2 2 2 2 
pfv  = 2y(d  -a  )u  + r(28+3y)(d  -a  )v 


(3.27) 


2 2 2 2 2 
pfu  = -d<j  + 2y(d  -a  )v  + r(28+3y)(d  -a  )u  + 2r(g+y)d  u (3.28) 


0 = a<J>  + 2u(d  -a  )w  - 2r(8+T)adu 


(3.29) 


dv  = - aw 


(3.30) 


where  d = d/dx  . Elimination  of  w between  (3.29)  and  (3.30)  gives 
2 

* 2u  22 

<p  = (d  -a  ) dv  + 2r(8+y)du 

a 


Substitution  of  this  in  (3.28)  gives 
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2U  2 2 2 2 2 
pTu  =>  - — - (d  -a  ) v + T C26  + 3y)  (d  -a  )u 


(3.31) 


Equations  (3.27)  and  (3.51)  constitute  a pair  which  should  be  solved 
for  the  unknowns  u and  v.  Introducing  nondintensional  distance  and  wave- 


number through 


C = 26 


k = a26 


and  writing  D = d/di;,  equations  (3.27)  and  (3.31)  become 


2 2 2o62r  T 22 

(D  -k  )u  = v-  (26+3y)  (D  -k  )v 


2 2 2 

(D  -k  ) v = 


— u + (26+3y)(D2-k2)u 


Making  the  transformation 


2 2 

V - jpj  .rk._  v 

u 


the  pair  (3.33)  becomes 


2 4 2 2 


22  4o  X r k o62r2k2  2 z 

(D  -k  )u  = ^ V - (26  + 3y)(D  -k  )v 


(D2-k2)2v  = - u + (2S  -y)  (D2-  k2)u 

4p6 


which,  on  introduction  of  the  non-dimensional  numbers 


2 4 2 

x = = 46^ 

V 

_ 26  + 3y 

H - 2 — » 

4p6 


(3.32) 


(3.33) 


(3.34) 


(3.35) 
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can  be  written  as 

2 2 2 2 2 2 
(D  -k  )u  = Tk  v - Tk  H(D  -k  )v  (3.36) 

2 2 2 2 2 

(D  -k  ) V = -u  + H(D  -k  )u  (3.37) 

Measuring  x from  the  lower  wall,  the  boundary  conditions  are 

u = v = Dv  = 0 at  C = 0 and  1 (3.38) 

Equations  (3.36),  (3.57)  and  (3.38)  define  a characteristic  value 
problem  which  has  solutions  only  if  T varies  with  k in  a definite  manner 
(for  fixed  H) . The  critical  value  of  T will  correspond  to  the  minimum 
of  the  T-k  diagram. 


3.9  Solution  of  the  Three-Dimensional  Disturbance  Problem 

The  characteristic  value  problem  (3 . 36) - (3. 38)  will  now  be  solved 
according  to  the  procedure  of  Chandrasekhar  (1954);  see  also  the  book 
of  Chandrasekhar  (1961),  pp.  300-303.  Since  u is  required  to  vanish  at 
C = 0 and  1,  it  can  be  expanded  in  a sine  series  of  the  form 


u 


C sinmirs 
m 


(3.39) 


Having  chosen  u in  this  manner,  we  next  solve  the  equation  obtained  by 
substitution  into  (3.37): 


2 2 

(D  -k  ) v = - 


C sinmir^ 
m 


+ 


2 2 
H(D  -k  ) 


C sinning 
m 


(3.40) 


m= 1 m= 1 

and  arrange  that  the  solution  satisfies  the  four  remaining  boundary 
conditions  on  v.  With  v determined  in  this  fashion  and  u given  by  the 
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assumed  series,  equation  (3.36)  will  lead  to  the  secular  equation 
for  T. 

The  solution  of  (3.40)  with  the  boundary  conditions  v = Dv  = 0 
at  ; * 0 and  1 is 


w 

= \ , m ?(H+  — — \ o/B^sinhks+A^ccoshkC+B^ssinhkc-sinmTTsi 

£—>  mV+k  mV+k2  k 1 2 2 J 


m=l 

where 


(3.41) 


BW 

2 


B(m)  _ _ m [k  + (_1)>nsinh|.j 

^ [sinh  k + (-l)mk  sinh  k] 


[sinh  k cosh  k - k + (-l)m(k  cosh  k - sinh  k) ] 
2 2 

A = sinh  k - k 


(3.42) 


Substituting  u and  v from  (3.39)  and  (3.41)  into  (3.36)  we  get 


. £ 


2 2 ; 

+k  )C  sinmirc*  Tk 
J m 


ao 

/ 2 2, 

1 — > m it  +k 


2 <H  + tK]  \ -2HkB(m)coshk; 
m tt  +k  l 2 


m=l  . 


m=l 


♦ (B^m^-2HkA^m^)sinhkC  + A^  Scoshk;;  + B ^ csinhkS  - [1+H(m  tt  +k  )]sinmTT(;> 

1 2 2 2 J 


Multiplying  the  above  by  sinm:rC  and  integrating  over  c = 0 to  1 results  in 

2 2 2 5 


Tk 


VA 
/ 2 2 

t — i m n + 


„ .(m  tt  +k  ) — 77  (H+ 

k2  ] 2Tk2  m2 


4-2^  -TTTI  {-2HkB 

tt  +k  ,n  tt  +k  '■ 


(m) 

2 


m=l 


[(-!)n+1coshkT-lJ  + (B(m)-2hkA(m))  (-l)n+1sinhk+(-l)n+1A(m)[coshk 


.2  2 

k +n 


— sinhk]+B^  [(-l)n+1sinhk - f (-l)n+1coshk+l ) 

ir  2 k +n  V 1 1 J 


- [1  + H(mV  ♦ k2)]  > = 0 


where  Smn  is  the  Kronecker  delta.  In  the  above  equation  the  requirement 

2 2 2 

that  not  all  of  the  coefficients  C / (m  it  +k  ) are  zero  is  the  vanishing 


of  the  following  infinite  determinent 


(m2ir2+k2)-^  + (H+  - j ■ ) -ny  ■ { -2HkB(m)[(-l)n+1coshk  + 1] 

| 2Tk2  mV+k2  nV+k  1 2 

+ (B^-2HkA^m^)  (-l)n+1sinhk+(-l)n+1A^m^  [coshk  - - ^ 2 i s^-n^  k] 

12  2 k +n  it 

+ B^m)[(-l)n+1sinhk-  ^Tz\ C-Dn+1sinhk+l)]|  - [l+H(mV+k2)  ]-^  J II  =0 


Substituting  the  values  of  the  constants  B|m^,  , B^  from  (3.42), 

the  above  simplifies  to 


CH+  1 i 2kmnn 

2 2 2'  222  2 2 
m n +k  (n  n +k  ) (sinh  k-k  ) 


2 2 2^  { ^Si 
n it  +k  '• 


inhkcoshk-k) 


\ 2 2 2 

X [l+(-l)n+mJ  + (sinhk-kcoshk)  [(-l)n+1+(-l)m+  ]>  +[1+H(m  ir  +k  )]  -HL 


2222  mn 

Mn  ir  +k  ) -20.  =0 

k T 


(3.43) 
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This  is  the  secular  equation  which  must  be  satisfied  in  order  that  our 
characteristic  value  problem  should  have  a solution.  A very  good 
approximation  to  (3.43)  is  obtained  by  setting  the  (1,1)  element  of  the 
matrix  equal  to  zero;  this  is  known  to  give  results  within  1%  accuracy 
for  Newtonian  fluids.  To  this  approximation  the  solution  of  (3.43)  is 


T = 


[l+H(ir2+k2)]2 


2 2 3 

(*  +k  ) 


16kTT2cosh2(k/2) 


(3.44) 


1 - 


(ir2+k2) 2 (sinhk+k) 


A plot  of  T vs.  k for  various  values  of  the  non -Newtonian  parameter  H 

will  now  be  obtained  from  (3.44).  To  choose  realistic  values  for  H = 
2 

(28+3y)/4p6  , note  that  8 and  y are  of  order  px  where  x is  some  time 
scale  of  the  fluid,  and  therefore 


H - rp  _ xT  v _ xf 
2p52  2 62r  2<5  + 2 

xT  being  the  ratio  of  the  time  scales  of  fluid  and  flow  is  taken  of 
order  unity,  and  6+  is  known  to  be  of  order  5.  Therefore  H - 1/50. 
Figure  3.5  shows  a plot  of  (3.44)  for  11=  ± 0.01.  From  the  stability 
diagram  it  is  seen  that  the  flow  is  stabilized  for  H < 0,  and 
destabilized  for  H > 0. 

It  may  be  noted  that  the  factor  involved  in  the  definition  of  H, 
namely  28  + 3y,  also  appears  in  Section  2.9  to  control  the  sign  of  the 
Weissenberg  effect:  H > 0 gives  positive  (free  surface)  V.’eissenberg 

effect  and  also  a destabilization  in  the  present  case. 

Note  that  the  curve  for  H = 0 (Newtonian)  is  valid  even  for  dis- 
turbances of  large  magnitude  since  the  nonlinear  terms  in  (3.24),  which 
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Figure  3.5  Neutral  stability  curves  for  disturbances  invariant  in  the 
flow  direction. 
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were  subsequently  dropped,  were  involved  in  the  terms  containing  B and  y 
only.  As  a matter  of  fact  the  curve  H = 0 yields  the  solution  of  Joseph 
(1966)  which  in  terms  of  the  sublayer  thickness  reads  5 = (1715/4) 4 = 

4.54.  The  agreement  with  Joseph's  solution  also  shows  that  the  approxi- 
mation made  by  setting  the  (1,1)  element  of  equation  (3.43)  to  zero  is 
quite  accurate. 


I 

3.10  Stability  of  Plane  Couette  Flow  of  Walters  B Fluid 


We  have  so  far  assumed  the  fluid  to  obey  the  second-order  constitu- 
tive equation.  For  two-dimensional  disturbances  we  have  seen  that  terms 
in  B vanished  and  the  results  indicated  stabilization.  Now  the  second- 
order  fluid  with  the  B-term  removed  is  identical  to  the  Walters  fluid 
A . However,  the  Walters  A fluid  predicts  a negative  Weissenberg  effect 
(when  the  free  surface  is  considered) , and  it  would  be  interesting  to 
investigate  the  present  problem  assuming  the  fluid  to  obey  the  Walters 
fluid  B : 


a.  . 

ij 


De.  . 

2“eij  - 


ekjVi,k 


- e..V.  ,) 

lk  j ,kJ 


[K  > 0] 


(3.45) 


which  shows  a positive  Weissenberg  effect.  Repeating  the  entire  pro- 
cedure, namely  that  of  maximizing  the  viscosity  for  all  disturbances 
which  have  instantaneously  stationary  total  energy,  we  found  that  the 
disturbances  obey  the  equation 


pE. .u.  = 
1J  J 


’V  [2y6ij+  K(-UkEij,k  + 4eikEk 


kj 


2Eikekj 


E,  .u.  , 
kj  i ,k 


Eikuj,k  + 6eikekj^ 


KE,  . ,-u. 
kj ,i  k, j 


(3.46) 


r 
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where  E. . and  e. . are  the  strain  rates  of  the  basic  and  disturbance 
ij  ij 

motions.  We  have  only  considered  the  special  case  of  plane  Couette 

flow  of  shear  rate  F in  presence  of  two-dimensional  disturbances. 

Defining  a stream  function  u = and  u = -ip  , elimination  of  $ 

1 ,2  2 ,1 

by  cross-differentiation  between  the  1 and  2 components  of  equation 
(3.46)  results  in 

4 

pr<p  + pV  + 2Kr^  = 0 (3.47) 

,12  >12jj 

which  is  identical  to  equation  (3.15)  derived  for  the  second-order  fluid. 
Therefore,  the  Walters  liquid  B also  predicts  a stabilization  of  a plane 
Couette  flow  with  two-dimensional  disturbances. 

3.11  Discussion 

The  fact  that  all  the  three  constitutive  equations,  namely  the 

t f 

second-order  fluid  as  well  as  the  Walters  fluid  A and  B , which  display 
different  types  of  normal  stresses  and  W'eissenberg  effects,  should  lead 
to  identical  results  in  two-dimensions  is  interesting.  Walters  (1962b) 
has  shown  that  the  two-dimensional  linear  stability  equations  for  the 
Walters  liquids  A*  and  b'  are  identical.  To  explain  all  this,  we  sub- 
stituted the  constitutive  equations  for  the  two  Walters  fluids  into  the 
equations  of  motion  for  the  total  (disturbed)  flow,  and  eliminated  the 
pressure.  In  both  cases  we  arrived  at  the  equation 

D 2 4 D 4 

pJ^Vip^uVip-K  — Vip  (3.48) 


1 

] 


1 


where  i (i  is  the  stream  function.  Moreover,  Chun  5 Schwarz  (196S)  have 
derived  equation  (3.48)  for  a second-order  fluid  in  two-dimensions,  in 
which  case  the  terms  involving  g dropped  out. 
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We  therefore  arrive  at  the  general  conclusion  that  in  two  dimensions 

» i 

the  second-order  fluid,  as  well  as  the  Walters  fluids  A and  B , all  give 
identical  results  concerning  the  velocity  distribution,  vorticity  dis- 
tribution, or  stability  characteristics;  their  normal  stress  behavior 
is,  however,  quite  different  from  each  other. 

The  assurance,  that  our  demonstration  of  stabilization  of  a plane 
Couette  flow  due  to  elasticity  is  not  merely  a result  of  a particular 
constitutive  equation,  is  gratifying. 


CHAPTER  4 


STABILITY  OF  PLAVE  POISEUILLE  FLOW  OF  OLDROYD  FLUID 


4.1  Introduction 


The  problem  of  the  stability  of  a plane  Poiseuille  flow  of  a 
Newtonian  fluid  has  been  solved  by,  among  others,  Nachtsheim  (1964). 

The  same  problem  has  also  been  solved  for  a second-order  fluid  by  Chun  § 
Schwarz  (1967),  for  a Maxwell  fluid  by  Tlapa  5 Bernstein  (1970)  and  for 
a Walters  fluid  by  Fong  5 Walters  (1965).  In  this  chapter  we  shall 
treat  this  problem  for  the  more  general  constitutive  equation 


+ X 


So 


St 


- 2y(e^  + X 


Sc 


St 


(4.1) 


which  has  been  called  "Oldroyd  Liquid  B"  in  Section  2.6.  Equation  (4.1) 
is  identical  to  equation  (2.18);  we  are  using  Greek  letters  for  the 
stresses  and  the  strain  rates  s1^  since  the  corresponding  Roman 
letters  are  being  preserved  for  the  basic  steady  motion.  The  Maxwell 
fluid  equation  is  obtained  from  (4.1)  by  setting  = 0. 

Unlike  the  case  of  plane  Couette  flow,  the  method  of  small  dis- 
turbances is  applicable  to  the  case  of  stability  of  a plane  Poiseuille 
flow.  We  shall  use  this  method  in  this  chapter. 


4.2  Nondimens ionnl  Equations  of  State  and  Motion 

In  the  rest  of  this  chapter  we  shall  use  only  cartesian  tensors. 

Using  the  rule  (2.14)^  for  expanding  the  convective  derivative  of  a 
doubly  contravariant  tensor,  the  equation  of  state  (4.1)  can  be  written  as 


- 


Ill 


a.  . + X (a..  . + V,  a..  , - a,  .V.  . - a.,V. 
ij  i ij,t  k ij,k  kj  i,k  ik  j,k' 


2ueij  * 2“Vsrj,t  * Vkei),k  - ‘kjVi,k  - EikVj,k) 


The  equations  of  continuity  and  motion  are  respectively 


V.  . = 0 

1,1 


p<Vi,t  + VkVi,k} 


= " w,i  + aik,k 


where  n is  the  pressure.  Let  26  denote  the  width  of  the  channel  and 
the  centerline  velocity  (Fig.  4.1).  Then  introducing  the  nondimensional 
variables 


X.  = 
1 

xt/6 

V.  = 
l 

V./U 
l m 

t = 

tUn/6 

IT  = 

TT/pii^ 

ij 

aij/pU: 

the  equations  of  state,  continuity  and  motion  become  respectively 


°ij  + A^aij,t  + Vkaij,k  " akjVi,k  “ aikVj,k^ 


R eij  + R ^Eij,t  + ^keij,k  Gkj^i,k  ~ eikVj,k^ 


(4.2) 


V.  . = 0 

1,1 


(4.3) 


V.  ♦ V,  V.  . = -it  . + a.  . 
i,t  k l , k ,i  ij,j 


(4.4) 
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where  the  tilde  (~)  on  all  symbols  have  been  omitted,  and  the  following 
nondimensional  numbers  have  been  introduced 


R = Reynolds  number  = 


pU  6 
m 


A = Relaxation  number  = 


X,U 
1 m 


(4.5) 


B = Retardation  number  = 


A2Um 


The  effect  of  nondimensionalization  was  to  eliminate  the  appearance  of 

-1 

p,  replace  X and  X by  A and  B respectively,  and  replace  y by  R 


4.3  Solution  of  the  Basic  Flow  • 

Let  the  disturbed  flow  be  written  as  a steady  basic  flow  plus  a 
time  dependent  disturbance,  assumed  small: 


vi 

= U.  + 

l 

u. 

i 

IT 

= P + 

P 

a.  . 
ij 

= s. . 

13 

+ s . . 
13 

e.  . 
il 

= E.  . 
13 

+ e.  . 
13 

(4.6) 


The  equations  of  state,  continuity  and  motion  for  the  basic  motion  are 

? 2R 

sij  ♦ ^Vij.rtA.k-^kV3  = r V ^ (ukEij.k  “ Vi.rtkV* 

U.  . = 0 

1.1 

0 = -P  . + S.  . 

.1  13.3 
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Assuming  a parallel  basic  flow  of  the  form  [U^ (x^) ,0,0] , a solution  of 
the  equation  of  state  yields  the  following  stress  matrix 


1,2 


u 


[S.  . ] = 

1 ir 


1,2 


Introduction  of  this  stress-tensor  into  the  equation  of  motion  gives 
the  result  that  P = constant,  and  that  the  velocity  profile  is 
parabolic: 

U = 2x  - x 2 0 < x <2  (4.7) 

1 2 2 ~ 2 ~ 


4.4  Stress  Change  Due  to  Disturbance 

Substituting  = S^  + s„  into  the  constitutive  equation  (4.2) 
for  the  total  motion,  subtracting  the  constitutive  equation  of  the  basic 
motion,  and  neglecting  nonlinear  terms,  we  get 


Sij+^Sij  , t+uk^i  j ,k+lJksij  ,k  skj^i,k~^kjUi,k  sikUj,k  Sikuj,k^ 


(4.8) 


2e.  . 


R R ^eij ,t+uk^ij ,k+^keij ,k  ekj^i,k  ^kjUi,k  eik^j,k  ^ikuj,k^ 


We  shall  assume  here  that  Squire's  theorem,  namely  that  two-dimensional 

disturbances  are  more  dangerous  than  three-dimensional  ones,  is  true 

for  this  problem;  Tlapa  5 Bernstein  (1970)  have  proved  it  for  plane 

Poiseuille  flow  of  a Maxwell  fluid.  Assuming  therefore  only  two- 

dimensional  disturbances  [u  =0,  ( ) =0],  the  significant  components  of 

3 » 3 


(4.8)  are 
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s + A(s  + U s ) 

11  11, t 1 11,1 


-lii  - Au  S + 2AU  s + 2A(S  u + S u ) 

R 2 11,2  1,2  12  11  1,1  12  1,2 


+ £“[u  + U u -(u  + u )U  -U  u ] 

R 1 l,lt  1 1,11  1,2  2,1  1,2  1,2  1,2 


s ♦ A(s  + U s ) 

12  12,  t 1 12,1 


u + u 

1,2  2,1 


- Au  S + AU  s + AS  u 

2 12,2  1,2  22  11  2,1 


(4.9) 


+ [%  (u  + u )+^uU  +%U(u  ♦ u ) 

R L V l,2t  2, It  2 1,22  1 1,21  2,11 


U U ] 
1,2  2,2 


s + A(s  + U s ) = 2AS  u + — K-2 

22  22, t 1 22,1  12  2,1  K 


+ — [u  +Uu  -U  u ] 

R 2,2t  1 2,21  1,2  2,1 


Now  assume  that  the  disturbances  are  periodic  in  the  x^-direction : 


u.  = u.  fx  ) e 
1 / 2 


ikx^ikct 


s.  . = s. . (x  ) e 

ij  lj"  2 


ikXj-ikct 


(4.10) 


p = p (x  ) e 


ikx^-ikct 


where  k is  a real  constant  — the  (nondimens ional)  wavenumber  of  the 
disturbance,  c = cr  + ic^  is  the  complex  wavespeed,  and  quantities  with 
(A)  are  complex  amplitudes;  as  usual  only  the  real  part  of  the  right 
hand  side  of  equation  (4.10)  represents  the  physical  quantity.  The 
motion  is  stable  or  unstable  according  as  c.  < 0 or  c^  > 0. 
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Substituting  (4.10)  in  (4.9),  and  defining 

3 (x  ) = 1 + ikA(U  -c)  (4.11) 

2 1 

we  find  that  the  stress  amplitudes  are  related  by 


2iku  „ - _ ~ . 

Bs  = n - - Au  S + 2AU  s + 2A(ikS  u + S Du  ) 

P5U  R 2 11,2  1,2  12  11  1 12  1 


X l*\  (Ux-c)  + Ui  2(2Dui  + iku2)] 


6s 


12 


DV  ikU2  . AU£_Ul^i  + au  s + ikAS  G (4.12) 
R R 1,2  22  11  2 


♦ X Vi  ik(Du^+  ikG2)  (U^c)  + % u1>22\-  U1>2Da21 


26 


2Du„ 


Us 


22 


,2  - —2  2B  r; 

ir  “2  ♦ 


IT  * X llkDWc)  - ikui>2u2] 


where  D = d/dx  is  a differential  operator. 

2 

/ 

4. 5 Equations  of  Motion 

Substituting  IT  = + u^^  into  the  continuity  equation  = 0, 

and  subtracting  from  it  the  continuity  equation  for  the  basic  flow 
Ui  i = 0,  we  get 

u.  . = 0 (4.13) 

i.i 

for  the  continuity  equation  of  the  disturbance.  Similarly  substituting 
(4.6)  into  the  equation  of  motion  (4.4),  and  subtracting  the  equation  of 
motion  of  the  basic  flow  from  it,  we  get 


■ 
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u.  . + u,U.  ,+U.  u.  ,=-p.  + s..  . 
- ,t  k i,k  k i,k  ij ,j 


(4.14) 


as  the  equation  of  motion  of  the  disturbance.  Introducing  periodicity 
(4.10)  into  the  continuity  (4.13)  and  the  two  equations  of  motion 
(4.14),  we  get 

iku  + Du  = 0 

1 2 


iku  (U  -c)  + u U = -ikp  + iks  + Ds 
11  2 12  r 11  12 


(4.15) 


iku  (U  -c)  = -Dp  + iks  + Ds 

21  12  22 


Substituting  the  stress  amplitudes  from  (4.12)  into  the  above,  the  two 
equations  of  motion  finally  become,  after  some  tedious  algebra  which 
results  from  the  coupling  of  the  stress  components  in  equation  (4.12), 

ikR8 [-  (U  -c)Dv  + vU  ] 

1 1,2 

2 3 2 

= k 3Rp  - { y (D  v-k  Dv)  + (26  -y  )D2v  + (6  + y )Dv 

,2  ,2  ,22  ,22 

2 8 . 2 2 

♦ 26  (6  -y  ) Dv}  + k y v + (y(D  v - k v)  + 2(6  - y )Dv 

*2  ,2  ,2  ,2  P ,2  ,2 


26  y 
, 22 ' 


+ (6  - y ) v + 2v6  (6  - y )}  + a 

,22  ,22  ,2  ,2  ,2  » 


Dv 


23  22(g  2'Y  2) 

■ v(6 - y ) ♦ ’-2  v 


,222  ,222 


6 


(4.16) 


„ 2 2 

ik6R(U  -c)v  = -8RDp  + y(D  v-k  v)  + 2(8  -y  )Dv 

1 ,2  ,2 


♦ (6  -y  )v  + 2v6  (8  - y ) 


,22  ,22 


,2  ,22  ,22 


(4.17) 
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where  we  have  written  u for  u , v for  u , and  defined 

1 2 


y(x  ) = 1 + ikB(U  -c) 
2 1 


(4.18) 


4.6  Generalized  Orr-Sommerfeld  Equation 


We  now  have  to  eliminate  p between  (4.16)  and  (4.17).  To  do  this, 

divide  (4.16)  by  S(x2),  differentiate  with  respect  to  x , multiply  it 

2 

by  8 again,  and  add  it  to  k times  (4.17).  This  yields,  after  some 
lengthy  algebra,  a generalized  Orr-Sommerfeld  equation 


2 2 

ikR8[(U  -c) (D  v-k  v)  - vU 

1 1,22 

4 3 2 

= vDv+bDv+bDv+bDv+bv 
3 2 10 


(4.19) 


where  the  coefficients  are  given  by 


\ = 2(5  ,(1  - V 

3 2 P 


(4.20) 


2 28 

b = -2yk  ♦ 38  (1  - J)  + 28  (8  -y  ) (1-  ^ (1  - J) 

2 ,22  P ,2, 2, 2 P P P 


(4.21) 


b = 28  (1  - J)  - 2k%  (1  - J)  + 48  (1  - i)(B  - y ) 

1 ,222  P ,2  P ,2  P ,22  ,22 


^8  2^  22  y ^2  1 

f1  ' }>  - -T-  Cl  - })  Wf2-  y>2) 


(4.22) 


b = yk  + 8 - y - k 8 (1-  i)  - 2k  8 (8  - Y ) 


* I K V * Q 

,2222  ,2222  ,22  P 


.2  ,2  ,2 


2 

2k  6 , 48  „ 36  „„ 

(y  8 " “a1  ‘ Y > g~  (8  - Y ) 

P ,2  ,2  P P ,222  ,222  P ,22  ,22 


(8  - Y ) 

6 2 »2  .2 


P 2 1 

(l  - £)  (8  - y ) 

a e ,22  ,22 


(4.23) 
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In  the  limiting  case  of  the  Newtonian  fluid,  we  have  X^  = X^,  which 
means  A = B and  therefore  3 = y;  the  four  coefficients  then  simplify  to 


b = 
3 

b = 
2 

b = 
1 


-23k 


3k 


and  (4.19)  reduces  to  the  Newtonian  Orr-Sommerfeld  equation 

ikR[(U  -c)  (0%  - k%)  - vU  ] = D*V  - 2k2Dv  + k4v 

1 1,22 

However,  on  substitution  of  y = 1 (Maxell  fluid),  our  coefficients 
turn  out  to  be  slightly  different  from  those  of  Tlapa  5 Bernstein  (1970); 
we  however  think  that  there  is  no  algebraic  mistake  in  our  calculations. 


4. 7 Simplification  for  Small  Relaxation  Times 

It  has  been  verified  by  Tlapa  5 Bernstein  (1968,  1970)  that  for 

values  of  the  stress  relaxation  number  A = X U /<5  in  the  range  of  0.0  to 

1 m 

0.1,  the  stability  diagram  is  virtually  unchanged  if  we  neglect  terms  of 
2 

order  A and  higher  in  the  generalized  Orr-Sommerfeld  equation.  This 
way  they  achieved  considerably  greater  speed  in  computation  without 
sacrificing  any  accuracy. 

It  should  be  mentioned  here  that  Tlapa  5 Bernstein  (196S,  1970) 

have  computed,  on  the  basis  of  the  experimental  data  of  Virk  et  al  (1967), 

-2 

that  A = X U /5  is  of  order  10  at  the  onset  of  the  Toms  phenomenon. 

1 m 

But  an  examination  of  the  procedure  of  Tlapa  5 Bernstein  (1968)  revealed 
that  it  was  erroneous;  what  they  computed  was  X ll^/D  where  D is  the 


A 
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diameter  of  the  pipe  and  is  the  average  velocity  in  turbulent  flow. 

However,  D/U^  in  turbulent  flow  has  the  significance  of  the  time  scale 

of  the  large  eddies,  which  cannot  have  any  fixed  relation  with  the  time 

scale  X of  the  molecules. 

1 

The  proper  time  scale  here  is  that  of  the  "wall  flow"  (all 

evidence  indicates  that  the  interaction  producing  the  Toms  phenomenon 

2 

occurs  in  the  wall  region),  which  is  v/u^  , where  v is  the  kinematic 

2 

viscosity  and  u^  is  the  friction  velocity  given  by  u^  = (wall  shear 

stress)/p.  For  the  two  fluids  that  Tlapa  f)  Bernstein  consider,  Virk 

2 

et  al  (1967)  give  A^u^.  /v  = 0.256  and  0.733.  In  considering  the 
stability  of  a laminar  motion,  the  time  scale  ratio  should  also  be  based 
on  the  wall  shear  stress.  In  a plane  Poiseuille  flow  this  gives 

-2 

2X  U /6  = 0.256  and  0.733.  Thus,  A is  considerably  larger  than  10  at 

1 m j 

the  onset  of  the  Toms  phenomenon. 

In  any  case,  in  the  present  work  we  shall  confine  ourselves  to 

values  of  A and  B in  the  range  0.0  to  0.1,  so  that  terms  of  order  equal 
2 2 

or  higher  than  A , B , AB  are  negligible.  Under  this  assumption  the 
coefficients  given  by  (4 . 20) - (4 . 23)  simplify  to 


b = 0 

3 

b = -2yk' 
2 

b - 0 

1 


(4.24) 


b = yk 
0 

and  therefore  (4.19)  reduces  to 

2 2 4 2 4 

ikRZ[(U  -c) (D  v-k  v)-U  v]  = D v - 2k  Dv  + k v 

1 1.22 


(4.25) 
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| 

I 


where 


Z(x  } 
2 


g 1+ikACUj-c) 

Y l+ikB(U^-c) 


(4.26) 


Equation  (4.25)  obviously  reduces  to  the  Newtonian  Orr-Sommerfeld 

equation  if  Z = 1,  and  to  the  simplified  Maxwellian  Orr-Sommerfeld 

equation  of  Tlapa  5 Bernstein  if  Z = 3. 

Equation  (4.25)  is  true  for  the  stability  of  all  two-dimensional 

parallel  flows  of  the  Oldroyd  fluid.  For  the  particular  case  of 

2 

Poiseuille  flow  the  basic  flow  is  given  by  U = 2x  - x 

-1  2 2 

The  boundary  conditions  are  that  the  disturbance  velocity 

components  u^  = u^  = 0 at  the  boundaries,  which  imply  v = Dv  = 0 at 

x = 0 and  x = 2.  However,  since  the  velocity  profile  is  an  even 

function  of  x^  about  the  line  x^  = 1,  the  disturbance  can  be  separated 

into  even  and  odd  function  parts.  The  former,  which  has  a simpler  flow 

pattern,  usually  gives  a lower  critical  Reynolds  number;  hence  the 

second  boundary  condition  at  x = 2 is  replaced  by  a condition  at  x =1 
3 

that  Dv  = D v = 0.  The  four  boundary  conditions  are  therefore 

v = Dv  = 0 at  x = 0 

3 2 (4-27) 

Dv  = D v = 0 at  x = 1 

2 

Equation  (4.25)  with  the  boundary  conditions  (4.27)  define  an  eigen- 
value problem;  for  given  k,  R,  A and  B a solution  (eigenfunction)  is 
possible  only  for  proper  values  (eigenvalues)  of  c.  For  constant  A and 
B,  it  is  of  interest  to  determine  the  minimum  of  R in  a R-k  diagram  on 
which  c^  = 0,  which  defines  the  margin  between  stability  and  instability; 


th  is  is  called  the  critical  Reynolds  number. 
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The  eigenvalue  problem  will  now  be  solved  by  the  numerical 
technique  developed  by  Nachtsheim  (1964),  who  solved  the  same  problem 
for  a Newtonian  fluid. 


4. 8 Initial  Value  Technique  of  Nachtsheim 

The  approach  to  the  eigenvalue  problem  for  fixed  k,  R,  A and  B is 
to  find  values  of  c = c^  + ic^  for  which  (4.25)  has  solutions  that 
satisfy  the  boundary  conditions. 

Trial  solutions  are  obtained  by  step-by-step  numerical  integration 
of  the  differential  equation  for  the  assamed  initial  values  and  an 
assumed  value  of  c.  The  proper  initial  values  and  c are  determined  by 
an  iterative  process  that  selects  the  one  solution  that  satisfies  the 
boundary  conditions. 

The  solution  is  started  at  x^  = 0 with  the  proper  boundary  values 

and  then  integrated  forward.  Another  solution  is  started  at  x =1 

2 

with  the  proper  boundary  values  and  then  integrated  backward.  Next  it 

is  necessary  to  perform  a matching  in  the  middle,  say  at  x£  = 0.5. 

For  computational  purposes  the  solution  is  carried  out  in  terms  of 

the  disturbance  vorticity  amplitude  and  the  x -velocity  amplitude.  A 

2 

measure  of  the  vorticity  amplitude  is  the  quantity  s defined  by  s = 

2 2 

(D  -k  )v.  To  show  this,  note  that  the  vorticity  is  given  by 


^2 

3x 


3u 

3x 


. ikx.-ikct  . 

JT  tv(Ve  ] ' 3T  [u(X2)e 

1 2 


ikXj-ikct 


= (ikv-Du)e 


ikXj-ikct 


h> 


The  vorticity  ainplitude  is  therefore  really  u = ikv  - Du;  the  quality 

s is  related  to  u through 

2 2 

s = D v-k  v = ik(-Du  + ikv]  = iku 

where  the  continuity  equation  iku  + Dv  = 0 has  been  used. 

Instead  of  solving  the  fourth  order  equation  (4.25],  we  can  there- 
fore solve  a system  of  two  second-order  equations 


■ i 2 

v = s + k v 

s"  = k2s  + ikRZ[(U^-c)s  + 2v] 


(4.28) 


where  instead  of  D we  are  now  using  prime  (')  to  denote  differentiation 

with  x : we  have  also  used  the  condition  D = -2  for  Poiseuille  flow. 

2 1,22 

The  boundary  conditions  are  v=v'  = 0 at  = 0,  and  v*  = s'  = 0 at 
x^  = 1.  For  details  of  the  numerical  procedure,  see  Nachtsheim  (1964). 


4.9  Equations  in  Real  Form 


The  real  differential  equations  are  obtained  by  separating  the 
original  equations  into  real  and  imaginary  parts.  For  example,  (4.28) 
written  in  real  form  are  as  follows: 


v"  = s + k v 
r r r 


n 2 

v.  = s.  + k v. 
ill 


(4.29) 


s = k s - kR{Z  f(U  -c  )s.-c.s  +2v.]  + Z.[(U  -c  )s  +c.s.+2v  1} 
r r rlt  i r'  l l r iJ  i1^  l r r l i rJ 


s.  = k s.  + kR{Z  [ (U  -c  )s  +c.s.+2v  ] - Z.[(U  -c  )s.-c.s  +2v.]} 
l l rL  i ry  r l l rJ  il  i rJ  l l r iJ 


The  differential  equations  for  the  variation  with  respect  to  the 
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eigenvalue  cr  for  both  the  forward  and  backward  solutions  are  needed  in 
this  method.  Differentiating  equations  (4.29),  these  are 


.<  2 

v = s + k v 
r,cr  r,Cr  r,cx 

••  2 

v.  = s . + k v. 

i,c  i,c  i,c 


= k s -aR(Z  [(U  -c  )s.  -c.s  +2v.  „ -s . ] 
r,cr  r,cr  rL ^ i tj  i,cr  l r,cr  i,cr 

+ Z [ (U  - c )s.-c.s  +2v. ] 
r,crL v i rJ  l i r iJ 

+ Z . [ (U  -c  )s  +c.s.  +2v  -s  ] 

ilv  i r,cr  l i,cr  r,cr  r J 

+ Z.  [ (U  -c  )s  +c.s.+2v  1) 
i,c  L i rJ  r l i rJ 

’ r 1 

ti  2 

s.  = k s.  +aR{Z  [(U  -c  )s  +c.s.  + 2v  -s  ] 

i,cr  i,cr  r 1 ^ i xJ  r,cr  i i,cr  r,cr  rJ 

+ Z f (U  -c  ) s +c . s . +2v  ] 
r,c  1 i r r l i rJ 

- Z. [ (U  -c  )s.  -c.s  + 2v.  -s.] 

i 1 r i»cr  1 r>cr  i»cr  i 

- Z.  [ (U  -c  )s.-c.s  +2v . ] } 

i,c  1 i ri  ir 


(4.30) 


In  equations  (4.29)  and  (4.30),  the  real  and  imaginary  parts  and 

their  derivatives  with  c of  the  non-Newtonian  function  Z(x  ) are 

r 2 


Z = 

r 


2 2 

(1+kAc.) (l+kBci)+k  AB(U  -c  ) 
(l+kBc.)2+k2B2(U  -c^2 


L 
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k(A-B) CU.-c  ) 

z.  = — — £ 

1 (l+kBc.)2+k2B2(U  -c^2 

2 

-2k  B(A-B) (U  -c  ) (1  + kBCj) 
’Cr  [(l+kBc.)2+k2B2(Ui-cr)2]2 


2 2 2 2 
(1+kBc . ) -k  B (U.-cJ 

Z.  = -k(A-B)  i — £ 

1,Cr  [ (1+kBc . ) 2+k2B2 (U  -c  )2]2 


(4.31) 


4. 10  Results  and  Comparisons 

The  computer  program  developed  by  Nachtsheim  for  the  Newtonian  case 
was  modified  to  take  into  account  the  presence  of  the  non-Newtonian 

parameters  A and  B.  The  program  has  been  stored  in  the  computer  library 

* 

of  The  Pennsylvania  State  University. 

The  computation  was  limited  only  to  the  region  near  the  minimum 

critical  numbers  R . The  forward  solutions  (started  at  x =0)  were 
cr  2 

matched  with  the  backward  solutions  (started  at  x = 1)  at  x =0.5. 

2 2 

The  range  x^  = 0 to  1 was  divided  into  128  steps  for  integration. 
Iterations  were  stopped  when  all  values  for  two  consecutive  iterations 
agreed  to  four  decimal  places. 

The  Newtonian  case  (A  = B)  was  programmed  first  to  see  whether  the 
results  agreed  with  those  of  Nachtsheim.  The  eigenvalues  calculated  at 
k = 1.00,  1.01,  1.02,  1.03  and  1.04  for  R = 5780,  5770  and  5760  agreed 


Off  campus  inquiries  requesting  source  decks  or  other  information  may 
be  directed  to:  Program  Librarian,  Computer  Building,  The  Pennsylvania 

State  University,  University  Park,  Pennsylvania  16802. 
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almost  exactly  with  those  of  Nachtsheim.  Fig.  4.2a  shows  plotted 
against  k for  various  values  of  R;  Fig.  4.2b  gives  versus  R 
corresponding  to  the  peaks  of  Fig.  4.2a.  Interpolation  gives  a minimum 
critical  Reynolds  number  of  Rcj>  = 5767  at  k = 1.02,  which  agrees  exactly 
with  Nachtsheim's  result. 

Eigenvalues  were  also  computed  for  five  non-Newtonian  liquids 
having  (A, B)  of  (0.02,  0.01),  (0.02,  0.0),  (0.01,  0.0),  (0.2,  0.1)  and 
(0.1,  0.0).  Figures  4. 3-4.7  show  the  determination  of  the  critical 
Reynolds  number  from  these  eigenvalues. 

The  critical  Reynolds  numbers  were  also  computed  for  the  same 
fluids  by  drawing  the  curves  of  neutral  stability  (c^  =0).  To  draw 
these  curves  (Fig.  4.8),  the  values  of  c^  were  written  in  a diagram  of 
R vs.  k,  and  points  of  c.  = 0 were  found  by  interpolation.  From 
Fig.  4.8  the  following  conclusions  can  be  drawn: 

(i)  The  presence  of  viscoelasticity  destabilizes  the  flow. 

(ii)  Comparing  the  curves  for  (0.02,  0.01)  with  (0.02,  0.0),  it 
seems  that  while  the  effect  of  stress  relaxation  alone  is  to  destabilize 
the  flow  compared  to  a Newtonian  fluid,  the  effect  of  strain  rate 
relaxation  at  a fixed  stress  relaxation  time  is  to  stabilize  the  flow. 
This,  of  course,  may  be  expected  since  the  Newtonian  results  should  be 
approached  as  A^  ■*  A^.  However,  it  should  be  realized  that  the  effect 
of  strain  rate  relaxation  alone  cannot  be  compared  to  a Newtonian  fluid 
since  the  inequality  A^  > A must  always  be  satisfied. 

(iii)  The  curves  of  (0.02,  0.01)  and  (0.01,  0.0)  coincide  almost 
exactly,  so  that  at  these  low  values  of  A and  B it  is  only  (A-B)  which 
is  important,  not  A and  B separately.  This  fact  could  also  be  brought 
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out  more  explicitly  if,  consistent  with  our  earlier  linearization  in 
A and  B,  we  had  also  linearized  equation  (4.26)  to  read  2=1+  ikX 
(A-B)  (Urc). 

(iv)  The  critical  wavenumber  increases  (that  is,  the  waves  become 
shorter)  as  the  fluid  elasticity  is  increased. 

The  conclusions  are  therefore  same  as  that  of  the  instability  of 
the  circular  Couette  flow  of  an  Oldroyd  fluid,  treated  in  the  next 

chapter. 

It  should  be  mentioned  here  that  our  results  for  the  Maxwell  fluid 
do  not  agree  with  those  of  Tlapa  5 Bernstein  (1968,  1970).  For  example, 
their  critical  Reynolds  numbers  for  A = 0.01  (B  = 0)  and  A = 0.1  (B  = 0) 
are  reported  [Tlapa  5 Bernstein  ( 196S) ] as  5767  and  5670  respectively, 
whereas  our  results  for  the  same  fluids  are  5745  and  5528  respectively. 
An  examination  of  their  computer  program  [Tlapa  5 Bernstein  (1968)] 
revealed  that  there  are  at  least  two  errors  in  their  program:  (i)  the 

presence  of  an  extra  factor  of  2 in  their  Runge-Kutta  integration 
formula  [the  22nd  card  of  their  subroutine  INTEGR  version  2 should  read 
CALL  F (X+H,  P,  R,  C3)  instead  of  CALL  F(X+H/2.,  P,  R,  C5) ] , and  (ii)  the 
absence  of  an  exponent  2 in  the  denominator  of  their  definition  of 
variable  DGR  (the  25th  card  of  their  subroutine  EQUTNS) . These  errors 
are  also  reflected  in  the  small  but  significant  inconsistencies  between 
Fig.  2 of  Tlapa  5 Bernstein  (1968)  [their  Newtonian  case]  and  Fig.  1 of 
Nachtsheim  (1964).  On  the  other  hand  our  Fig.  4.2  agrees  exactly  with 
Nachtsheim's  Fig.  1. 

In  addition  to  the  above  errors  in  Tlapa  8 Bernstein's  work,  there 
are  other  inconsistencies  between  our  algebra  and  theirs,  as  already 
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mentioned  in  Section  4.6.  However,  these  inconsistencies  vanish  for 
the  simplified  case  of  small  A and  B,  for  which  their  results  as  well 
as  ours  are  valid. 

4.11  Discussion 

From  the  present  as  well  as  the  other  studies  mentioned  in 
Section  1.3(vi),  it  appears  that  all  the  stability  investigations  on  a 
plane  Poiseuille  flow  have  predicted  a destabilization  due  to  elasticity, 
although  the  constitutive  equations  assumed  display  different  types  of 
normal  stress  behavior  and  the  Weissenberg  effect.  This  seems  to  be  a 
result  of  the  two-dimensionality  of  the  problem,  for  it  was  demonstrated 
in  Section  3.11  that  in  two  dimensions  the  vorticity  equations  for  the 
entire  (disturbed)  flow  are  identical  for  the  second-order  fluid  and 
the  Walters  short  memory  fluids  A*  and  b'.  Although  this  fact  cannot 
probably  be  demonstrated  in  general  for  the  Oldroyd  fluids  (for  which 
the  stresses  cannot  be  explicitly  expressed  in  terms  of  the  strain  rates), 
it  must  be  true  at  least  for  the  small  relaxation  times  assumed  in  the 
present  study;  this  is  because  at  small  relaxation  times  the  Oldroyd  A 
and  B fluids  con  be  approximated  by  the  Walters  short  memory  fluids  a' 

f 

and  B respectively,  as  shown  in  Section  2.7. 


5.1  Introduction 

The  stability  of  circular  Couette  flow  has  always  been  a happy 
hunting  ground  of  stability  investigators.  One  reason  is  that 
analytical  solutions  can  be  obtained  without  a great  deal  of  difficulty. 
Moreover  the  theoretical  results  can  be  easily  compared  with  experiments 
for  this  geometry.  From  the  discussion  of  Section  1.3  it  can  be  seen 
that  this  problem  has  not  been  treated  for  the  Oldroyd  fluid.  However 
we  think  that  the  Oldroyd  constitutive  equation  is  the  most  reliable  of 
all  the  simple  constitutive  equations  commonly  used.  For  this  reason 
the  stability  of  a circular  Couette  flow  of  the  Oldroyd  fluid  B is 
investigated  in  this  chapter.  The  small  disturbance  theory  is  applied. 

5 . 2 Some  Properties  of  Cylindrical  Coordinates  and  General  Tensors 

In  this  chapter  we  shall  use  the  notation  of  general  tensors,  with 
covariant  suffixes  written  below  and  contravariant  suffixes  above,  and 
a comma  shall  denote  a covariant  derivative  instead  of  simply  a space 
derivative . 

The  cylindrical  polar  coordinates  chosen  are  denoted  by 

1 2 3 

x=r,  x = 0,  x = z (5.1) 

whose  metric  tensors  are 


The  scale  factors  defined  by  = /g^  (no  sum)  are  therefore 


h = 1 

1 


h = r 

2 

h = 1 

3 


and  the  Christoffel  tensors  can  be  shown  to  be 

X 


r = -r 
22 

2 2 

r = r 

12  21 


r 


T*  = 0 otherwise 


(5.3) 


(5.4) 


The  physical  components  of  a vector  or  tensrr  will  be  noted  by  suffixes 
enclosed  within  parentheses.  The  rule  for  transforming  a mixed  second- 
order  tensor  b*  into  its  physical  components  b^j^  is  given  by 


(ij) 


h .b 
i J 


h. 

J 


(no  sum) 


(5.5) 


and  the  rules  for  transforming  other  types  of  tensors  can  be  guessed  from 
the  above.  Thus,  for  example,  b^  = h^b*  = b^/lw  for  3 vector.  The 
covariant  derivative  of  a mixed  second-order  tensor  is 
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b?r\ 

J Pk 


(5.6) 


and  the  expressions  for  covariant  derivatives  of  other  types  of  tensors 
can  be  deduced  from  this. 

A'i  the  above  rules  of  tensor  analysis  will  be  used  in  the  next 
few  sections.  For  a simple  treatment  of  these  rules  see,  for  example, 
the  book  of  Aris  (1962). 


5.3  Governing  Equations 


The  equations  of  continuity,  motion  and  state  are  respectively 


V1.  = 0 

,1 


(5.7) 


,3V  . li  ij 

p(-rr-  + v v •)  = -ir  -g  J + a J . 

>J  »J 

ij  , <5o^  - , ij  . <5e^, 

a J + X -t— — = 2p(e  J + X -j- — ) 

l 6t  v 2 ot 


(5.8) 


(5.9) 


Here  V1  is  the  velocity,  e.  . = %CV.  . + V.  •)  is  the  strain  rate,  a1-*  is 

7 ij  i,J  J,i 

the  extra  stress,  m is  the  pressure,  p is  the  density,  \i  is  the  viscosity, 

\ is  the  stress  relaxation  time  and  X is  the  strain  rate  relaxation 
1 2 

time  (also  called  'retardation  time').  The  notation  6/5t  denotes  the 
convective  derivative  and  is  given  by  the  rules  (2.14),  that  is, 

k 

operating  on  a mixed  tensor  b^,  it  is  given  by 

6bk  3b^  . . 

* Vmbk  + Vm.bk  - Vk  bm  (5.10) 

6t  3t  i,m  ,i  m ,m  i 


Substituting  the  rule  (5.6)  for  the  covariant  derivative  into  (5.10),  it 
can  be  easily  shown  that  the  terms  involving  the  Christoffel  tensors 
vanish,  giving  simply 
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U = 0 

r 


U0  - rU 


uz=  0 


ur  = o 


u9  = ” 

r 


u = o 


The  covariant  components  of  the  strain  rate  tensors  E^ 
are  then  determined  as 


(5 .13) 


[E.  .] 


- — A 

2 dr  V 


1 A A 

2 dr  V 

0 


from  which  the  contravariant  components  E1^  = g^g^E 


g g'  are  found  to  be 


[Elj] 


1 A A 

2 dr  ^rJ 


1 A A 

2 dr  V 


(5.14) 


Substitution  of  the  above  into  the  constitutive  equation  of  the  basic 
flow 


cij  . .aU1  ckj  3UJ  cik.  , _ij  , , ,3UX  _kj  3UJ  rik 

S - A ( S J + — - S ) = 2uE  2uA  ( E J+  E 


1 3xk 


3x" 


yields  the  stress  components 

[sij]  = 


0 

d ,U, 


2 3xk- 


d ,U, 

M dF  (F> 


ax'1 


U (-)  2p  (A  ~A  ) [A  (~)]2 
dr  r 12  dr  r 1 


0 

0 


(5.15) 


141 


Substitution  of  the  above  stress  matrix  into  the  equation  of  motion 

0 = - P .g1;i+  S1^  . 

>J  ,1 


yields  the  velocity  distribution 


where 


U = Er  + — 
r 


R2n2-Rin, 

r2-r2 
2 1 

RKcni-°2) 

2 2 

R -R 
2 1 


(5.16) 


(5.17) 


Here  and  R are  the  radii  of  the  inner  and  outer  cylinders  respec- 
tively, and  and  are  their  angular  velocities  (Fig.  5.1). 


5.5  Stress  Change  Due  to  Disturbance 


Let  us  confine  our  discussion  to  the  state  of  neutral  stability  so 
that  3/3t  = 0,  and  let  the  physical  components  of  the  velocity  of  the 


disturbed  motion  be 


hr)’  “ 

he)  - “ * » 

hr)  ’ “ 


(5.18) 


where  u,  v and  w are  functions  of  r and  z.  Substituting  o1-1  = S1-1  + 
s1-*  into  the  constitutive  equation  (5.S)'  for  the  total  motion,  sub- 


tracting the  constitutive  equation  of  the  basic  motion,  and  neglecting 


the  nonlinear  terms  we  get 


3u  u 3w 

— + — + — — = 0 (5.21) 

3r  r 3z 


Subtracting  the  equation  of  motion  for  the  basic  flow  from  that  of  the 
disturbed  total  flow,  neglecting  the  nonlinear  terms  and  setting  3/3t=0 
for  neutral  stability  we  get  the  equation  of  motion  of  the  disturbance: 


-1_ 


A 
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p(U3uX  . + uJUL  •)  = - p .g13  + s 13 

>J  >J  >J 


(5.22) 


The  three  components  of  the  above  equation  give 


♦ -4v  - 

r 

11  11 
9p  9s  s 

9r  + 9r  r 

13 

9s  22 

+ - ~ rs 

9z 

(5.23) 

2pEu 

r 

12 

9s 

9r 

12 

3s 

+ 

r 

23 

9s 

9z 

(5.24) 

0 = , 

9z 

13 

+ s'  + 

9r 

13 

+ 

r 

33 

9s 

9z 

(5.25) 

Using  the  expressions  (5.20)  for  the  disturbance  stresses,  the  above 
equations  respectively  reduce  to 


”2p  (E  + 

r 


u[ 


2 

9 u 
9r2 


+ 


I iE 

r 9r 


2 

u 9 u, 

T + 

r 9z 


+ 


8y(ArA2)F 


2 

r 


— (-) 
9r  V 


24F  UA1(A1-A2) 


5 

r 


9u 

9 


(5.23)' 


2pEu= 


1 9v 
r 9r 


2u(ArA2)F 

_ 

r 


r9  u 

[_T  + 

9r 


i iy 

r 9r 


0 


I ly. 

r 9r 


+ 


9z 


(5.24) 


(5.25)  ' 


To  simplify  the  analysis  we  shall  make  the  usual  assumption  that  the 
gap  6 between  the  cylinders  is  small  compared  with  their  radii.  IVe  now 
introduce  the  nondimensional  radial  coordinate  measured  from  the  center 


of  the  channel  (see  Fig.  5.1) 
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r-R, 


C = 


C5.26) 


and  the  kinematic  parameter 


1 


(5.27) 


Assuming  that  ^ is  small,  it  can  be  shown  from  the  expression  (5.17) 
for  the  constants  in  the  basic  velocity  distribution  that 


E = 


Rl^l 

26 


F = - 


Riafi1 

"26 


E + = a [1  ♦ a(C  + %)] 

2 1 


Under  this  small  gap  approximation,  equations  (5.21),  (5.23)  -(5.25) 
become 

(5.28) 


1 3u  3w  _ . 

3z 


6 3C 


2 2 

, . . , 13p  ,1  3 u 3 u., 

-2pfi^  [l+a(i;+J2)  ] v-  - g-  + u(—  — 2 + 2 “ 


1 32u  32u,  4U(VW  H 

,2  3; 


,2  „ 2 2 

6 3£  3z 


6pA1  (X1-Xz)Sl1  Rxa  5u 

,3  H 


(5.29) 


pR  S^au 


2 

, 1 3 v 

u(  2 2 

6 3C 


0 = 


3 v 


) + 


u(A1-A2)Riaia 


3z 


I?  * -4 

o 


2 

3 w 
3^2 


2 2 
.1  3 u 3 u. 

(“2  7T  + 7T) 

6 3£  3z 


32w, 

3z2 


(5.30) 


(5.31) 
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5.7  Stability  Equations 


Now  assume  that  the  disturbances  are  periodic  in  the  z -direction 
with  wavenumber  a 

u = u(Osinaz 
v = v(c)sinaz 

(5.32) 

w = w(C)cosaz 
p = p(^)sinaz 


Introducing  the  nondimensional  wavenumber 

k = a6 

the  equation  of  continuity  (5.28)  gives 

w = t-  Du 
k 

and  equation  (5.31)  gives  the  pressure 

u 2 2 

p = D(D  -k  )u 

5k 


(5.33) 


(5.34) 


(5.35) 


where  D = d/dc.  The  remaining  equations  (5.29)  and  (5.30)  then  become 
2 2 

2pfl  J 6 k 2 2 2 2 

[1+a (£+%)]v  = (D  -k  ) u + 4 (A  -A  )Q  ak  Dv 

V 12  1 


6A  (A  -A  )n  2R  aV 

1 1 2 1 i 

+ ; DU 


(5.36) 


pR  flau  = T (D  -k  )V  + — l (D  -k  )u  (5.37) 

11  6 


Making  the  transformation 

2 2 

2pQ  5 k 

o 1 o 

v -+  v 

V 


equations  (5.36)  and  (5.37)  finally  become 
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222,  , 2 , 

[l+a(C+'i)]v  = (D  -k  ) u+2 (A-B) aDv-3A(A-B) ak  TDu 


(D2-k2)v  = -Tk2[u  - (A-B)  ([)2-k2)u] 


(5.38) 

(5.39) 


where  the  following  nondimens ional  numbers  have  been  defined 


T = Taylor  number  = - 


2 3 2 

2aG  6 RlP 


A = Nondimens ional  relaxation  time  = 


B = Nondimens ional  retardation  time 


yX 

p6 

uX 

p6 


The  boundary  conditions  are 


v = u = Du  = 0 


at  e = ±% 


(5.40) 


(5.41) 


Equations  (5.38)  and  (5.39)  subject  to  (5.41)  constitute  a character- 
istic value  problem;  that  is,  a nontrivial  solution  is  possible  only  if 
T varies  with  k in  a special  manner  (for  given  A,  B,  and  a). 

Equations  (5.38)  and  (5.39)  were  also  derived  by  Thomas  5 Walters 
(1963)  for  the  case  of  Walters  fluid  b',  given  by  equation  (2.25),  of 
which  the  Oldroyd  fluid  B is  a special  case.  However,  the  general 
solution  was  not  given. 


5 . 8 Solution  of  the  Characteristic  Value  Problem 

The  characteristic  value  problem  will  now  be  solved  as  follows: 
assume  first  an  expansion  for  u in  terms  of  a set  of  orthogonal  functions 


which  automatically  satisfy  the  four  boundary  conditions  u = Du  = 0 at 
£ = ±%.  This  expansion  is  then  substituted  into  (5.39),  which  then  can 
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be  solved  for  v,  subject  to  the  remaining  two  boundary  conditions  on  v. 
These  results  are  finally  substituted  into  (5.3S),  from  which  we  obtain 
that  a secular  determinant  of  infinite  order  must  vanish  in  order  to 
have  nontrivial  solutions. 

The  crucial  point  in  this  method  lies  in  the  choice  of  the  initial 
approximation  for  u.  A good  choice  of  the  orthogonal  functions  used  in 
the  expansion  of  u will  mean  a smaller  size  of  the  secular  determinant 
that  must  be  taken  in  order  to  achieve  a desired  accuracy.  Chandrasekhar 
§ Reid  (1957)  [see  also  the  book  of  Chandrasekhar  (1961),  p.  634] 
suggested  the  use  of  the  orthogonal  functions 


cnCO  - 


snU) 


coshin; 

cosh'iA 


sinhu  c 
n 

s~inh’2p 


cosA  ? 
n 

cos'jA 


sinv 

sin-'su 


(5.42) 


where  A^  and  are  the  positive  roots  of  the  equations 


tanh  rr  + tan  % A = 0 


coth  H U - cot  H u = 0 


The  functions  have  the  orthogonality  properties 


,/  cmU)cn(c)dc  - J s„(t)Sn(0*it  ■ 


J Cn(C)Sn(t)d5  - 0 

my2 


(5.43) 


(5.44) 


They  form,  therefore,  a complete  set  of  orthogonal  functions  in  the 


interval  (--’i,  H)  • 
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These  functions  have  been  used  by  Reid  (195S)  for  solving  the 
stability  of  a Newtonian  fluid  flowing  through  a curved  channel  due  to 
a pressure  gradient;  Thomas  $ Walters  (1963)  extended  that  work  to  a 
Maxwell  fluid.  Thomas  5 Walters  (1964)  also  applied  these  functions 
for  solving  the  instability  of  circular  Couette  flow  of  a Maxwell  fluid. 
We  shall  apply  the  method  here  to  the  Oldroyd  fluid.  It  may  be 
mentioned  here  that  since  Thomas  5 Walters  (1964)  omitted  all  the 
analytical  steps  and  reported  only  the  resulting  graphs,  we  shall  not 
have  a check  until  we  arrive  at  the  final  numerical  results. 

We  therefore  assume  a general  expansion 


[A  c (C)  t B s (;)] 
L m m m m J 


(5.45) 


Substituting  this  in  (5.39)  we  get 


2 2^  2 
(D  -k  ) v = -Tk 


) {A  [C  - (A-B) (C"  - k2C  )] 
£ | 1 mL  m v 1 K m m J 


It  2 

+ B [S  - (A-B) (S  - k S )]} 
ml  m J K m nr  J J 


(5.46) 


which  has  to  be  solved  subject  to  v = 0 at  g = ± The  homogeneous 


solution  is  obviously 


vu  = K coshkg  + K sinhkC 
H 1 2 


If  a particular  solution  is  assumed  as 


\ , * Ml  I II  III 

/ IE  C +F  C +G  C +H  C +e  S +f  S +c  S +h  S t 

l_ , 'mmmmmm  m mmmm6mm  mmJ 
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then  substitution  into  the  differential  equation  and  comparison  of 


coefficients  gives  (k  / A^  or  u^) 


4 4 


Tk  A 2 

m [k  - (A-B)  (Am  - k )] 


A 4-k4 
m 


F = 0 


m 


2 

Tk  A 


A 4-k4 
m 


H = 0 


m 


m 


Tk  B 2 44 

m [k  - (A-B)  Cum  - k )] 


4 V4 

y _ -k 
m 


f = 0 


2 

Tk  B 


m 


4 u4 
"k 
m 


h = 0 


ra 


The  total  solution  is  therefore 


v = K coshk^+K  sinhk^-Tk 
1 2 


m=  1 


-rhr  fc,lk2-  t*-6’  <C-k') 

A -k 
m 


+ C } - Tk 
m J 


I 


B 2 4 4 i i 

— — {S  [k  - (A-B)  (u  - k )]  ■*■  S } 
4,4‘m1  m raJ 


m=l 


ym  -k 
m 


The  constants  K and  K are  now  determined  from  the  boundary  conditions 
1 2 


v = 0 at  ; = t'i.  The  solution  finally  becomes 
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v = Tk 


T — 1 A 2\m  2 4 4 " 

) !2—  { £r-r  coshk;-C  [k  -(A-B)(X  -k  )]-C  } 

/ .4,4  1 cosh  jk  mL  mrn 

L i X -k 

m=l  m (5.47) 

CO  2 

’ B 2y„  2 4 4 >t 


- £ 


2 — {-.m-  sinhk^-S  [k  -(A-B)(y  -k  )]-Sj 

4 , 4 1 sinh'2k  mL  mm 

y -k 
m 


Substituting  the  expressions  for  u and  v given  by  (5.45)  and  (5.47)  into 
equation  (5.3S)  we  get 


(1+ja+aOTk 


EA  2X  2 4 4 It 

— { rn-  coshkc-C  [k  -(A-B)(X  -k  )]-C  I 

4 4 1 cosh^k  mL  mm 

A “ K 

m=l  m 


liSk  - s«» 

y_  -k 
m 


SACA-BJaTkV^.B^) 


A 2X  k 
m r m 


4 4 


. 2(A-B)aTk  — sinhkc-^ 


-Cjk  - (A-B)  (1  - k )1  - C ) 


X 4-kn 
m 


B 2y  k 
m i m 

+ 4,4  ' sinh’jk 


2y  k i 2 44  iii. 

2 — coshk;  - S Jk  - (A-B)  (u  -k  ) ] - S } 


y_  -k 
m 


(5.48) 


Now  multiply  the  above  equation  by  Cn(0  anc*  Sn(0  in  turn,  integrate 
between  g = -H  to  h,  and  use  the  orthogonality  properties  (5.44).  The 


following  two  algebraic  equations  result: 


L . 

m=l 


r ~ v."  - >1 


, 4 4 1 cosh  2k 

X -k 
m 


m J mn  mn 


1,44  2 

- _i_  rx  + k ) 6 - 2k  X J 

9 l <■  m J mn  mn 


Tk 


E 


m=l 


2u 


4 ,4  lsinh:ik 
Um  'k 
m 


(rSr  « - [E-  (a-bjo/V))  J01  - J2it 


4 4 


* 3A(A-B)aI  - t^W  P ' ^ ' W'8’  ' k’)I1„-  ho1 

1 U u -k 
m 


= ( 


and 


(5.49; 


w 

E 


m=l 


2X  2 4 4 

± { 2—  V - [k  - (A-B)  (X  -k  )]I  -I  } + 3A(A-B)aJ 

4 , 4 lcosh>2k  L m 01  21 


X -k 
m 


2X  k 2 4 4 

IVzllZ  { jl—  Q - [k  - (A-B)  (X  - k )]  J - j } 

L.  Li  InmchliW  x 1 k vm  10  30 


,4,4  lcosh'2k 
X -k 
m 


E 


m=l 


1 + ; m 


2y_  2 44 


4 ,4  ‘sinhjk 

_k 
m 


irrsrar  0 - lk  - - k ^Ein  - ¥mn> 


1 , 4 4 2 

-±-  {(U  - k )6  - 2k  Y } 

2 1 V-Mm  mn  mnJ 

Tk 


= 0 


(5.5C 


where  we  have  defined  the  following  functions  of  k: 


~ — — 


153 


P(k) 

= / coshkcC  ds 

i n 

Q(k) 

= / sinhk^S^dc 

U(k) 

J/2 

h 

= / ssinhkcC  d; 

V (k) 

. — *5 

H 

= / scoshk;S  d; 

following 

— 2 

constants  independent 

of  k 

I 

= / cc  S dc 

J 

% 

= / cs  Cd? 

01 

, m n 

-'■i 

01 

-Jj  m n 

I 

= / ;C  S dc 

J 

* " 

■ , / tsmcndt 

21 

t.  ® n 

21 

I 

= / SC  d; 

J 

' 

= / C S dc 

10 

v m n 

10 

i * m n 

I 

= / sm  Cd; 

J 

• 

= / C S dt 

30 

t/  ® n 

“"Z 

30 

, ^ m n 

X 

" 

= / C C d; 

Y 

> ■■ 

■ J S»Sndt 

-Jj 

mn 

k m n 

“"2 

mn 

(Note  that  the  symbols  for  these  integrals  are  not  really  meaningful  or 
systematic;  some  of  them  correspond  to  those  used  by  earlier  workers.) 
The  integrals  P(k),  Q(k),  U(k)  and  V(k)  can  be  evaluated  explicitly  with 
the  following  results: 


P(k) 


1 


cosh 


sinh';(>.  -k)  sinh'iCX  +k) 

•f 


X -k 
n 


X -t-k 

n 


COM  A 


ksinli’^kcos’.X  +X  coshhksinhX 

n n n 


, 2 .2 
X + k 
n 


(5.51) 
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T 


Equations  (5.49)  and  (5.50)  can  be  written  in  the  abbreviated  form 


\ {a  [a  -JLB]  + bc}  = o 

/ 1 ml  mn  Tk.2  mnJ  m mn 

m=l 
00 

]T  < Van  * BJE™  - ^2  Fml  1 ■ 0 


where  we  have  defined  the  following  functions  of  k: 


(5.49) 


(5.50) 


• 7T7  Sifer  P - t'2-  - xmn 

A "K 

m 

4 4 2 

B = (A  + k )6  „ - 2k  X 
mn  m mn  mn 


2u  2 4 4 

: = ■ t~  — . *T  U - [k  - (A-B)  (u  -k  )]J  - J 

mn  4 , 4 sinh-ik  m 01  21 

um  -k  V 

+ 3A(A-B)  a I - ILV-II*  f P-[L2-(A-B)(p^4  )]1  -I 

Blai10  4 4 I sinh-ik  L m 10  30 

Um  \ 


2A  2 44  \ 

D = - r TT  v - [k  - (A-B)  (A  -k  ) ] I - I 

umn  4 k4  cosh^k  L m 01  21 I 

m " \ ' 

2 

* “CA-BJal,-  ^Q-tk2-CA-B)Cx;.k4)]Jio-J 


(5.55) 


1°  A H-k“ 

m 


Wm  _k 
m 


10  30 


4 4 

F = (U  + k )<5  - 2kY 

mn  m mn  mn 
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V 


It  is  apparent  that  as  n varies  from  1 to  equations  (5.49)'  and 

I 

(5.50)  represent  a doubly  infinite  system  of  linear  homogeneous 

equations  for  the  coefficients  A and  B . The  condition  that  the  con- 

stants  A and  B do  not  all  vanish  is  that  the  determinant  of  the  system 
mm  1 

should  vanish,  that  is 


A 


mn 


1 


D 

mn 


B 

mn 


C 

mn 


E 

mn 


0 (5.56) 


This  is  the  required  characteristic  equation  from  which  the  curve  of 
neutral  stability  can  be  obtained:  Equation  (5.56)  yields  multiple 

roots  of  T,  but  only  the  smallest  positive  root  of  T need  be  considered. 


5.9  Numerical  Calculations 


The  simplest  approximation  to  the  characteristic  equation  (5.56)  is 
obviously  a second  order  determinant,  which  corresponds  to  taking  m = 
n = 1,  that  is  when  one  uses  one  even  and  one  odd  function  in  the  initial 
representation  for  u.  For  Newtonian  fluids  this  approximation  is  known 
to  give  results  within  1%  of  the  exact  values,  and  we  shall  assume  that 
the  accuracy  of  the  method  will  remain  good  for  our  viscoelastic  case 
since  we  shall  confine  ourselves  with  small  values  of  the  viscoelastic 
parameters  A and  B. 

We  shall  therefore  approximate  (5.56)  by 


A - — r B 

11  TV2  11 


'll  Tk2  11 


4 2 2 

(A  E -C  D )k  T - (B  E +A  F )k  T + B F =0  (5.57) 

11  11  11  11  11  11  11  11  11  11 


The  two  roots  of  T as  a function  of  k were  calculated  from  the  above 
equation,  and  the  smaller  positive  root  of  T(k)  defines  the  neutral 
stability  curve. 

The  values  of  the  various  constants  for  m = n = 1 are 


X 

= 

4.730041  p 

= 

7.853205 

1 

1 

I 

- 

0.1478245  J 

- 

0.1478245 

01 

01 

I 

- 

-0.00060475  J 

= 

-6.682854 

21 

21 

I 

- 

3.34202  J 

- 

-3.34202 

10 

10 

I 

= 

-122.06503  J 

— 

122.06503 

30 

30 

X 

- 

-12.30262  Y 

= 

-46.05012 

11 

11 

Some  of 

these  values 

were  already  evaluated 

by 

Reid  (1957)  and  Thomas  5 

Walters 

(1963).  The 

ne’.v 

■ ones  were  evaluated  by 

direct  integration  and 

then  checked  by  numerical  evaluation  of  the  integrals  in  the  computer. 


5.10  Results  and  Conclusions 


For  computing  numerical  results  two  cases  were  considered: 

(i)  Cylinders  rotating  in  the  same  direction  with  the  same  angular 
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velocity,  i.e.,  a = 0,  and  (ii)  Outer  cylinder  stationary,  i.e.,  a = -1. 

The  values  of  the  viscoelastic  parameters  A and  B were  varied  from  0.0 

to  0.1.  A few  sample  neutral  stability  curves  are  displayed  in 

Figs.  5.2  and  5.3.  All  the  results  are  finally  shown  in  compact  form 

in  Figs.  5.4  and  5.5,  which  plot  the  critical  Taylor  number  T£  against 

A,  with  B regarded  as  a parameter.  Figures  5.6  and  5.7  likewise  display 

the  variation  of  the  critical  wavenumber  k as  a function  of  A and  B. 

c 

The  results  for  B = 0 agree  with  those  of  Thomas  5 Walters  (1963)  for 
the  Maxwell  fluid. 

From  the  above-mentioned  figures  the  following  observations  can  be 

made : 

(i)  The  presence  of  viscoelasticity  destabilizes  the  flow. 

I 

(ii)  At  a fixed  X , an  increase  of  X (<  X ) stabilizes  the  flow. 

1 2 1 

This  of  course  must  happen  since  the  results  should  approach  their 

Newtonian  values  as  X -*■  X . 

2 1 

(iii)  From  Fig.  5.4  (a  = 0)  it  is  seen  that  each  curve  for  B = 

constant  is  exactly  identical  to  the  others,  only  displaced  horizontally. 

(Numerical  values  show  it  more  explicitly,  e.g.,  T for  A = 0.02,  B = 

C 

0.01  is  exactly  the  same  as  the  T for  A = 0.03,  B = 0.02  - both  being 
equal  to  1426.14.)  This  means  that  it  is  only  (A-B)  which  counts,  not 
A and  B separately.  That  this  must  happen  can  be  seen  from  an  examina- 
tion of  the  stability  equations  (5.38)  and  (5.39),  which  for  a = 0 show 
that  only  the  viscoelastic  terms  in  (A-B)  remain  in  the  equations. 

From  Fig.  5.5  (a  = -1)  it  can  be  seen  that  although  the  curves  do 
not  exactly  coincide  when  displaced  horizontally,  they  are  however 
close  to  being  so.  (For  example  T^  for  A = 0.02,  B = 0.01  is  2821.44, 


NONDIMENSIONAL  WAVE  NUMBER  k 


Figure  S. 


Neutral  stability  curves  for  some  pairs  of  values  of 
(A, B) , with  a = *0.0. 
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I 

while  Tc  for  A = 0.05,  0.02  is  2517.99.)  This  means  that,  even  for 
a t 0,  for  all  practical  purposes  it  is  only  (A-B)  which  is  important, 
not  A and  B separately.  However,  this  fact  is  not  expected  to  be  true 
for  higher  values  of  A and  B when  the  last  term  on  the  right  side  of 
equation  (5.38)  will  become  larger. 

(iv)  From  Figs.  5.6  and  5.7  it  is  apparent  that  the  presence  of 
viscoelasticity  increases  the  critical  wavenumber,  that  is  decreases 
the  cell  size,  at  the  onset  of  instability. 

All  the  above  theoretical  conclusions  are  therefore  identical  to 
those  of  the  plane  Poiseuille  flow  of  an  Oldroyd  fluid,  studied  in  the 
last  chapter. 

A discussion  of  the  experimental  observations  of  the  Couette  flow- 
instability  was  given  in  Section  1.4,  from  which  it  would  be  seen  that 
the  predictions  of  the  destabilization  and  the  increase  of  critical 
wavenumber  due  to  viscoelasticity  agree  with  at  least  half  the  observers. 


CHAPTER  6 


SOME  UNSTEADY  PARALLEL  FLOWS  OF  OLDROYD  FLUID 

6 . 1 Introduction 

This  chapter  is  not  directly  related  to  the  rest  of  the  disserta- 
tion. For  the  sake  of  curiosity  some  unsteady  parallel  flows  of  the 
Oldroyd  fluid  were  solved,  and  these  will  be  recorded  in  this  chapter. 
Sometimes  no  great  effort  will  be  made  to  interpret  the  results. 

Some  unsteady  parallel  flows  of  the  second-order  fluid  have  been 
solved  by  Ting  (1963)  and  Markovitz  fi  Coleman  (1964);  Tlapa  5 Bernstein 
(1968)  have  solved  the  problem  of  the  suddenly  accelerated  plate  (Stokes 
first  problem)  in  a Maxwell  fluid;  however,  no  such  works  seem  to  have 
been  done  for  Oldroyd  fluids. 

6.2  Governing  Equation 

The  differential  equation  governing  the  velocity  distribution  in 
unsteady  parallel  flows  of  Oldroyd's  liquid  B will  be  derived  in  this 
section.  The  velocity  components  in  a parallel  flow  are 


U = U (x  ,t) 

1 1 2 

U = U = 0 

2 3 


(6.1) 


The  constitutive  equation  for  Oldroyd  liquid  B was  given  in 


equation  (2.18);  in  cartesian  coordinates  it  is 
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S..  + A (S . . „ * U.  S . . . -S..U.  , - S . . U , ) 
ij  1 ij.t  k ij,k  kj  1 , k ik  j ,k^ 


2pE.  + 2uA  (E.  . + U.  E.  . - E,  .U.  , - E..U.  . ) 

ij  2 1J>t  k ij,k  kj  i,k  ik  j,k' 


(6.2) 


where  S.  . is  the  extra-stress,  and  E.  . = h([).  . + U.  .)  is  the  rate  of 

iJ  ij  v i,J  J,iJ 

strain.  For  the  case  of  parallel  flow  (6.1),  the  various  components  of 
the  constitutive  equation  (6.2)  give 


2 


+ Ms,, 

4 

■ U s ) 

= 2A 

s u 

- 2uA  (U  ) 

(6. 

,3a) 

1 11 

1 11,1 

1 12  l,: 

2 2 1,2 

+ \ (S 

-♦ 

■US  ) 

= UU 

+ yX  U + A S U 

(6. 

3b) 

1 12 

1 12,1 

1,2 

2 1 , 2t  1 22  1,2 

S 

+ X 

(S  + U 

S 

) = X 

S U 

(6. 

,3c) 

13 

1 

13,  t 

1 13 

,1 

1 23  1,2 

s 

+ A (S 

+ U 

s ) 

= 0 

(6. 

,3d) 

22 

1 22, t 

1 22,1 

s 

+ X (S 

+ u 

s ) 

= 0 

(6. 

,3e) 

23 

1 23, t 

1 23,1 

S 

+ A (S 

+ u 

s ) 

= 0 

(6. 

,3f) 

33 

1 33, t 

1 33,1 

Using  the  transformation 

£ = x -U  (x  )t 

1 1 2 

h = t 


equations  (6.3d),  (6.3e)  and  (6.3f)  can  be  expressed  as 


S.  . + A S . . 
ij  1 ij.n 


(ij  = 22,  23,  33) 


whose  solution  is 

S.j  = AU)e"n/Xl  = A(x  -U  t)e“t/Xl 
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If  we  imagine  that  the  flow  started  from  rest  , then  A(x^  - 0 and, 


therefore. 


S = S = S =0 
22  23  33 


Equation  (6.3c)  then  shows  S = 0,  so  that  the  stress  tensor  has  the 


[S]  = S 


S S 0 

11  12 


0 0 


0 0 


From  physical  considerations,  the  derivatives  of  the  above  extra- 


stresses with  respect  to  are  zero. 

Now  examine  the  equation  of  motion 


'ui,« 


-P  . + S. . 

.1  1J,J 


(6.4) 


Consideration  of  the  x - and  x -components  of  this  eauation  shows  that 

2 3 

P cannot  vary  in  these  directions,  i.e.,  P = P(x  ,t).  Next  consider  the 


x -component 


-P  = pU  - S 
,1  l,t  12,2 


The  left  hand  side  can  be  a function  of  (x^,t)  only,  while  the  right 
hand  side  can  be  a function  of  (x  ,t)  only.  Therefore  we  must  have 


■P  = pU  - S = f(t) 
,1  l.t  12,2 


(6.5) 


Differentiating  equation  (6.3b)  with  respect  to  x2,  setting  ^ 


Alternatively,  etc.  being  proportional  to  expf-t/A^)  will  decay 
very  rapidly  because  ^ is  a very  short  time. 
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S = 0,  and  using  (6.5),  there  results  the  following  differential 
equation  governing  the  velocity  distribution: 


pU  + pX  U = pU  + pX  U + f + X f 

l,t  1 l,tt  1,22  2 l,22t  1 »t 

If  there  is  no  pressure  gradient  in  the  direction  of  flow, 
and  therefore 


(6.6) 

then  f = 0 


3U. 


2 

3 U, 


'1  , ’ “1 

P TT~  + pX  7 = p 

3t  1 3t2 


2 

3x  2 


3 

3 U, 


+ pX 


2 3x  3t 
2 


(6.7) 


If  X = X , then  the  above  reduces  to  the  Newtonian  diffusion  equation 

1 2 


V 


2 

3 U 


3x 


1 

2 


2 


where  v = p/p.  Equation  (6.7)  may  be  compared  with  the  corresponding 
equation  for  the  second-order  fluid 


3t 


= P 


2 

3 U, 


3x 


3 

3 U 


1 


3x  23t 
2 


derived  by  Ting  (1963),  and  Markovitz  5 Coleman  (1964).  Note  that  in 
both  these  cases  the  order  of  the  differential  equation  has  increased 
due  to  the  presence  of  viscoelasticity. 


6 . 3 Oscillation  of  One  Plate  Parallel  to  Another 

Consider  an  Oldroyd  fluid  undergoing  a nonsteady  shearing  flow 
between  two  parallel  plates  (Fig.  6.1),  one  of  which  is  at  rest  at 
y = 0,  while  the  other  located  at  y = 5 undergoes  sinusoidal  oscillations 
in  the  x-direction  with  velocity  UQcosu)t.  The  pressure  gradient  in  the 
direction  of  flow  will  be  assumed  zero.  This  problem  has  been  solved 
by  Markovitz  5 Coleman  (1964)  for  a second-order  fluid. 
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IVe  seek  a solution  of  equation  (6.7)  of  the  form 

iut 


U^y.t)  = U(y)e 


(6.8) 


where  we  are  using  y in  place  of  x^.  Actually  only  the  real  part  of  the 


right  hand  side  is  meant,  i.e.. 


U = ILcosujt  - UTsinu)t 

1 R I 


(6.9) 


where  the  complex  amplitude  is  U(y)  = + iU^.  The  boundary  conditions 

read 


U(0)  = 0 

0(6)  = U 


(6.10) 


Substituting  (6.8)  into  (6.7),  and  using  the  boundary  conditions  (6.10), 
the  solution  obtained  is 

0 . u (6.11) 

0 sinhkd 

where 


k = 


iwfl  + iuX^ 


v(l  + ia>A 


(6.12) 


Denoting  k = a + ib,  the  above  gives 

A 


a = 


b = 


I+cj^X  , 2X  2+u2  (X  2+X  2)  - <o(A  -AJ 
1 2 ^ 1 2 1 2J 


|2v(1+uj2X22)  ^ 


/WU^A^2^2^2)  + u(xrx2) 


(6.13) 


From  (6.12),  the  real  and  imaginary  parts  of  U become 

U. 


UD  = — [sinhay  cosby  sinha6  cosb5  + coshav  sinby  coshad  sinbol 

k q 


tL 


(6.14) 


Uj  = — [coshay  sinby  sinhad  cosbd  - sinhay  cosby  coshad  sinbd] 


where 
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2 2 2 2 

q = sinh  a5cos  b<5  + cosh  afisin  bo  (6.15) 

The  solution  is  therefore  (6.9),  along  with  (6 . 13) - (6 . 15) . For 
different  values  of  the  non-Newtonian  parameters  uA ^ and  uA  , and  of 
the  quantity  3/u>/2v  , the  solutions  at  various  nondimensional  times  wt 

are  shown  in  Figs.  (6. 2) -(6. 6).  The  quantity  5 /to/ 2 v can  be  rewritten 

2 l' 

as  [(6  /2v)/(l/w) ] 2,  and  therefore  has  the  interpretation  of  square  root 
of  the  ratio  (viscous  diffusion  time) /(time  scale  of  flow);  for  small 
values  of  this  quantity.  Fig.  6.4  shows  that  the  velocity  distribution 
is  linear. 

Compare  the  velocity  distributions  for  a fluid  with  small  elasticity 
[say  uA^  = 0.2,  wA^  = 0.1  (Fig.  6.3)]  and  a fluid  with  high  elasticity 
[wAj  = 2.0,  gjA^  = 0.2  (Fig.  6.6)].  In  the  former  case  the  velocity 
distributions  for  eitner  wt  = 0 or  tot  = m/2  cross  the  y-axis  only  once, 
whereas  in  the  latter  case  they  cross  the  y-axis  two  times;  the  highly 
elastic  case  is  therefore  beginning  to  resemble  a wave  propagation 
through  an  elastic  solid. 

By  solving  (6.3a)  and  (6.3b),  and  of  course  assuming  that  the  flow- 

started  from  rest,  it  can  be  shown  that  the  stress  components  are  given 

by  2iwt 

2h(^1-^2)e 

SH  = (1  + ituX  )(l  + 2iwA  ) 

4 * * 

1+iuA,  . . , , 

_ ...  2 itot  coshkv 

S = uU  k — : — — e ■■ 

12  0 1+iuA  smhka 

Thus,  the  shear  stress  oscillates  with  the  same  frequency  (w)  as  that  of 
the  velocity,  while  the  normal  stress  oscillates  with  2w;  this  fact  is 
also  true  for  a second-order  fluid  [Markovitt  S Coleman  (1964)]. 


UQkcoshky 

sinhk6 


1.0 
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Velocity  profiles  for  wX1  = 2.0,  w>  2 = 0 2,  ■- 
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Now  let  us  consider  the  special  case  of  d/u/2v  very  snail,  that  is, 

when  the  viscous  diffusion  time  is  much  smaller  than  the  time  scale  of 

the  flow;  this  case  will  correspond  to  any  of  the  following:  (i)  small 

gap,  (ii)  slow  oscillation,  or  (iii)  highly  viscous  liquid.  Under  this 

condition  ao , bd,  and  hence  k<$,  will  all  be  much  smaller  than  unity. 

Then  the  hyperbolic  sines  in  (6.11)  can  be  expanded  according  to  the 
3 5 

rule  sinhx  = x + x /3!  + x / 5 ! + ...,  and  keeping  only  one  term  of  the 
series  we  get  sinh  ky  : ky,  sinh  kd  = kd,  so  that 

0 . V 


Therefore,  from  (6.9),  we  get 


u(y.t) 


V 


coswt 


which  agrees  with  the  linear  profile  of  Fig.  6.4. 

The  shear  stress  in  this  limiting  case,  found  from  (6.36),  is 


12 


V (!♦<*>  \X2) 
6(l+u2*  ’) 


(coswt  + 4) 


where 


* -tan 


-1 


w(ArX2) 


1+u  \ \ 

1 2 


The  shear  stress  therefore  oscillates  with  the  same  frequency  as  that  of 
the  velocity,  but  is  out  of  phase  with  it  by  an  angle  4>  which  approaches 
zero  as  -*■  X . It  is  likely  that  this  pliase  difference  could  be 
related  to  the* stored  energy  of  the  fluid. 
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6.4  Oscillation  of  a Single  Plate  in  an  Infinite  Fluid 

Consider  a flat  plate  oscillating  harmonically  parallel  to  itself 
in  a semi-infinite  fluid;  this  is  sometimes  called  "Stokes'  Second 
Problem,"  and  its  Newtonian  solution  is  treated  in  standard  texts  [see, 
for  example,  Schlichting  (1968),  p.  85].  Iv'e  shall  however  consider  an 
Oldroyd  fluid,  for  which  the  governing  differential  equation  is  (6.7), 


subject  to 


U (o.t)  = U costot 

1 0 


and  also  that  U + 0 as  y -*■  ®.  Keeping  the  Newtonian  solution  in  mind, 
let  us  seek  a solution  of  the  form 


U^y.t)  = Uoe"kycos[wt-<Hy)]  (6-1 

where  <j>  is  the  phase  difference  at  height  y;  we  are  going  to  determine 
k and  <f>(y)  • Subst  tuting  (6.17)  into  (6.7),  we  get 


2 

d <f> 

-tO  -V  r- 

dy 

2 

♦ -A  to  - 
1 


+ 2vk  j vA  to  f— 7“}  + vA  cok  sin(oot-<J>) 

dy  2 ldy J 2 


vk  + v {4^-}  - vA  to  ^-1  + 2vA  <ok  cos(tot-<j>)  = 0 
dyJ  2 dy2  2 dy 


which  can  be  true  for  all  t only  if 

2 2 

-to  -v  — j + 2vk  ^ - vA  to  {^-}  + vA  tok  = 0 

a„z  dy  2 Ldy J 2 


2 2 2 2 

•A  to  -vk  + v {■— } - vA  to  ^-1  + 2vA  tok  = 0 

1 dy  2 d„2  2 dy 


Elimination  of  (d$/dy)  between  these  two  equations  gives 


11  - 2k  ^ = 
dy2  dy 


1+A  A to 
to  1 2 

V ,,22 

1+A  to 
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The  solution  of  this  equation  is 

2 

2 kv  1 ♦ A, Aw 

» „ *K/  w 1 2 

# = A + Be  + ttt—  — — - y 

T 2kw  , ,22  7 

1 ♦ A„  w 
2 

where  A and  B are  constants.  The  condition  $ = 0 at  y = 0 requires 
B = -A,  so  that  the  solution  is 


2ky  u 1 + 

♦ " - e )*  jfc  — “¥r 

1 ♦ A u 
2 


(6.21) 


To  determine  A,  we  require  that  A = should  result  in  the  Newtonian 
expression  for  <J>,  which  is  [see,  for  example,  Schlichting  (1968), 
p.  85)]  <pN  = wy/2vk.  This  requires  A = 0 in  equation  (6.21),  giving 


1 + A Aw 
u 12 

2kv  . ,22 

1.+  A w 
2 


(6.22) 


To  determine  k,  the  above  ..xpression  for  $ is  substituted  into 
(6.19),  whose  solution  which  agrees  with  the  Newtonian  solution  at 

xi  = Xz  is  §iven  by 


-w  (Aj-A2)  + w/w2(A1-A2)2  + (1  + AiA2w  )' 


2 2 

2v(l  + Az  w ) 


(6.23) 


The  positive  sign  outside  the  braces  ( } was  chosen  so  as  to  insure 

2 

Uj  / 00  at  y = 00 . Using  A > it  can  be  easily  shown  that  k >0. 

The  velocity  distribution  given  by  (6.17),  (6.22)  and  (6.23)  is  plotted 
in  Figs.  6. 7-6. 9 for  different  values  of  the  non-Newtonian  parameters 
and  for  various  nondimensional  times  wt.  (Note  that  our  Newtonian 
solution  does  not  agree  with  the  numerical  values  in  Schlichting's 
Fig.  6.8.  His  ordinates  seem  to  be  /T  times  their  correct  values.) 
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Figure  6.7  Velocity  distribution  in  the  neighborhood  of  an  oscillating 
wall.  Comparison  of  a Newtonian  fluid  with  one  having 
wA^  = 0.4,  wA^  = 0.1. 
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Compare  the  velocity  distributions  for  a fluid  with  small 

elasticity  [say  u>A  = 0.4,  toA  =0.1  (Fig.  6.7)]  and  a fluid  with  high 
^ 2 

elasticity  [toA^  = 2.0,  toA^  = 0.2  (Fig.  6.9)].  In  the  former  case  the 
velocity  distributions  for  either  tot  = 0 or  tot  = ir/2  cross  the  y-axis 
only  once,  whereas  in  the  latter  case  they  cross  the  y-axis  three  times; 
the  highly  elastic  case  is  therefore  beginning  to  resemble  the  wave 
propagation  through  an  elastic  solid. 


6.5  Shear  Waves  in  Oldroyd  Fluid 

It  is  well-known  that  in  incompressible  fluids  without  a free 

surface,  the  only  type  of  possible  wave  motion  is  that  of  the  so-called 

"shear  waves"  or  "transverse  waves"  which  propagate  perpendicular  to 

the  fluid  velocity.  To  see  this,  assume  a wave  motion  given  by  the 

velocity  vector  . , . 

7 l (a  x -ut) 

■ I n mm' 

U.  = C . e 

J 3 

where  denotes  the  components  of  the  wavenumber  vector.  Tne  divergence 
of  the  velocity  is  then 

i(a  x-wt) 

■■  . „ 'mm 

U.  . = lC.a.e  = lU.a. 

3>J  J J J 3 

The  condition  of  incompressibility  IL  ^ = 0 therefore  requires  that  the 
vectors  U and  a must  be  perpendicular  (Fig.  6.10). 

Now  assume  a shear  wave  in  an  Oldroyd  fluid,  of  the  form 


U 


. , iay-ito  t to.  t 

Aei(ay-»t)  , A „ t 'l 


(6.24) 


where  a is  real  but  « = u ♦ ito. . Substituting  (6.24)  into  the  governing 
equation  (6.7),  we  get 
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2 2 2 
iw  ♦ u A - va  + iva  A u>  = 0 
1 2 


(6.25) 


This  gives  ui  = 0 (and  hence  the  wave  speed  c = 0)  if  v = 0,  confirming 
the  fact  that  shear  waves  cannot  propagate  in  an  ideal  liquid.  The 
solution  of  (6.25)  is 


-i(l+va  X^)  + vAX 1 va2- (l+va^Ap ■ 
— 


. 2 \ >2 


CO  = 


(6.26) 


Two  cases  are  possible  depending  on  the  sign  of  the  quantity  under  the 
radical  sign. 

2 2 2 
Case  I:  4A  va  > (1  + X vc  ) 

1 2 

Then  (6.26)  shows  that  u>  > 0 and  io.  <0,  so  that  the  waves  are 

r l 

propagated  as  well  as  damped  in  time. 

2 2 2 
Case  II:  41  va  < (1  + X va  ) 

1 2 

Then  (6.26)  shows  that  u>r  = 0,  < 0,  so  that  the  waves  do  not  propagate 

but  are  "standing"  and  damped,  as  in  Newtonian  flows. 

2 

It  is  not  difficult  to  find  a meaning  of  the  inequality  4A  va  > or 

2 2 

< (1  + X2va  ) . We  have  come  across  this  inequality  in  many  other 

parallel  flows  of  the  Oldroyd  fluid  that  we  tried,  even  in  motions  not 

apparently  of  the  wave  propagation  type  (in  which  cases  a signified  the 

Fourier  components  into  which  a general  velocity  distribution  can  be 

expanded).  We  suggest  the  following  interpretation  of  the  inequality, 

2 2 2 

which  for  Case  II  is  easily  recast  in  the  form  4va  (X  -X  ) < ( 1 - X va  ) . 
First  note  that  the  Newtonian  fluid  satisfies  this  inequality.  As  A^  is 
gradually  increased  at  constant  A^,  this  inequality  is  satisfied  in  the 


m 


■ 
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beginning  but  not  for  large  X ^ X . Thus,  Case  II  represents  small 
elasticity,  whereas  Case  I represents  large  elasticity. 

We  can  therefore  conclude  that  for  small  elasticity  the  waves  are 
purely  diffusive  and  are  therefore  simply  damped  in  time.  On  the 
other  hand,  a fluid  having  large  elasticity  develops  some  "elastic 
solid"  like  characteristics,  so  that  the  waves  are  also  propagated. 

For  the  special  case  of  a Maxwell  fluid  (X  = 0),  the  above 
phenomenon  can  be  given  a more  complete  explanation  as  follows.  From 

the  model  of  the  Maxwell  fluid  of  Section  2.12  it  is  easily  seen  that 

★ 

the  shearing  modulus  of  the  fluid  is  G = 2y/X  ^ . The  speed  of  propaga- 
tion of  elastic  waves  is  then  CE  = /G/p  = /2v/X  . On  the  other  hand 
the  speed  of  viscous  diffusion  is  Cv  - va.  One  can  consider  that  the 

viscoelastic  disturbances  are  propagated  if  Cj;  > Cy,  and  are  simply 

2 2 

damped  if  Cj;  < Cy.  This  gives  XjVa  > 2 for  propagation  and  X^va  < 2 
for  simple  damping,  which  agree  with  the  inequalities  under  Cases  I and 
II,  except  for  a numerical  factor  which  could  be  avoided  by  defining 

Cv  = I'll  va. 


Lumley,  in  a private  communication,  has  shown  an  alternative  way  of 
finding  an  expression  for  it.  For  a rapidly  applied  stress  the 
Maxwell  fluid  displays  the  characteristics  of  a purely  elastic  solid, 
and  gives  X^  = 2ue,  where  S is  stress  and  e is  displacement.  Integra- 
tion gives  XjS  = 2ue,  from  which  G = 2u/X  . 


CHAPTER  7 


SUMMARY  AND  CONCLUSIONS 


Three  problems  on  the  stability  of  viscoelastic  flows  have  been 
solved.  In  addition,  some  consideration  is  given  to  the  energy 
significance  of  the  common  viscoelastic  constitutive  equations,  some 
interesting  observations  have  been  made  on  the  influence  of  the  method 
of  measurement  on  the  sign  of  the  Keissenberg  effect,  and  some  unsteady 
parallel  flow  problems  of  viscoelastic  fluids  have  been  solved.  A 
summary  of  the  results  and  conclusions  is  given  below. 


7.1  Plane  Couette  Flow 

It  has  been  argued  that  the  application  of  the  small-disturbance 
linear  theory  to  the  investigation  of  the  stability  of  a plane  Couette 
flow  of  a viscoelastic  fluid  is  incorrect.  Consequently  the  energy 
method,  following  the  classical  variational  technique  of  Orr,  has  been 
used  for  solving  the  plane  Couette  flow  problem  of  a second-order  fluid 
obeying  the  constitutive  equation 


(2) 


= yA.  . + BA. , A,  . + yA. 
lj  ik  kj  ' ij 


where  S. . is  the  extra-stress,  A- • is  twice  the  strain  rate  tensor, 

ij  ij 

(2) 

Ajj  is  the  second  Rivlin-Ericksen  tensor,  y is  the  viscosity  in  simple 
shear,  and  S and  y(<0)  are  two  elastic  properties  of  the  fluid. 

It  has  been  found  that  in  order  for  the  energy  method  to  become 
applicable  to  a viscoelastic  fluid,  a new  concept,  namely  that  of  the 
elastic  potential  energy  of  a viscoelastic  fluid,  has  to  be  introduced 
and  an  expression  for  it  has  to  be  discovered. 
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No  assumption  regarding  the  size  of  the  disturbances  was  made  for 

the  case  of  two-dimensional  disturbances,  for  which  the  nonlinear  terms 

vanished  from  the  stability  equations  and  a simple  analytical  solution 

could  be  obtained.  The  results  (Fig.  3.3)  indicate  that  the  presence  of 

viscoelasticity  stabilizes  a plane  Couette  flow  of  this  fluid.  For 

example,  for  a Newtonian  fluid  the  critical  Reynolds  number  is  found  to 
2 

be  Rcr  = pFS  /u  = 44.3,  where  ? is  the  shear  rate  of  the  basic  flow,  p 
is  the  density,  and  5 is  half  the  distance  between  the  plates.  On  the 
other  hand  for  a viscoelastic  fluid  with  G = YT/u  (which  is  really  the 
negative  of  the  ratio  of  the  time  scale  of  the  fluid  to  the  time  scale 
of  the  flow)  = -0.5,  it  is  found  that  Rcr  = 57.8. 

It  was  suggested  by  G.  I.  Taylor  that  half  the  distance  between  the 
plates  at  the  critical  state  in  a plane  Couette  flow  is  a good  measure 
of  the  thickness  of  the  viscous  sublayer  in  any  turbulent  flow.  Noting 
that  R = pT6  /y  = u^p  6 /y  = (6  ) , where  ut  = /ur/p  is  the  friction 
velocity,  we  get  the  result  that  the  nondimensional  viscous  sublayer 

thickness  in  a Newtonian  fluid  is  <5+  = /44 . 3 = 6.65,  whereas  for  a 

■ 

+ , 

viscoelastic  fluid  with  G = -0.5, it  is  5 = /57.8  = 7.6;  this  corresponds 

to  an  increase  of  14.3%  in  the  sublayer  thickness.  IVe  have  thus 
suggested  an  explanation  of  the  experimentally  observed  fact  that  when 
small  quantities  of  polymer  molecules  are  added  to  a turbulent  pipe  flow 
of  an  ordinary  Newtonian  liquid,  a thickening  of  the  viscous  sublayer 
occurs;  this  in  turn  results  in  an  increase  of  the  flow  rate  at  the 
same  wall  stress  (the  Toms  phenomenon). 

For  the  two-dimcns ’ onal  case  mentioned  above,  the  terms  involving 
8 had  dropped  out  of  the  final  equations,  and  therefore  the  results 


1 
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would  also  be  valid  for  the  Walters  fluid  a'  (which  is  in  fact  identical 
to  the  second-order  fluid  with  the  5-term  removed).  However,  we  also 
solved  the  two-dimensional  problem  for  the  Walters  fluid  b',  which 
displays  a completely  different  (and  more  realistic)  type  of  normal 
stress  behavior  and  the  Weissenberg  effect.  The  results  were  identical. 
The  prediction  of  identical  results  for  these  three  constitutive 
equations  has  been  shown  to  be  a consequence  of  two-dimensionality;  it 
has  been  demonstrated  that  the  two-dimensional  vorticity  equation  for 
the  entire  (disturbed)  flow,  which  is  obtained  after  elimination  of 
pressure,  is  identical  for  all  the  three  fluids. 

The  same  problem  was  also  tried  with  more  realistic  type  of  dis- 
turbances, namely  those  in  the  form  of  longitudinal  rolls  invariant  in 
the  flow  direction.  Only  the  second-order  fluid  was  investigated  for 
this  type  of  disturbances.  The  resulting  equations,  however,  involved  a 
large  number  of  nonlinear  terms,  and  the  problem  was  linearized  by 
assuming  that  the  disturbances  were  small.  The  validity  of  this  process 
of  linearization  in  an  application  of  the  energy  method  is  questionable; 
however,  it  should  be  noted  that  the  linearization  was  necessary  only  in 
the  elastic  terms,  and  the  solution  is  probably  not  meaningless.  The 
results  (Fig.  3.5)  indicate  that  the  critical  Reynolds  number  may 
increase  or  decrease  due  to  the  presence  of  elasticity  in  the  fluid, 

2 

depending  on  the  sign  of  the  nondimensional  quantity  H = (2S  + 3y)/4p<5  . 

It  may  be  noted  that  the  case  H > 0 which  gives  a destabilization  in 
this  flow  also  gives  a positive  Weissenberg  effect  (observed  by  the  free 
surface  in  a circular  Couctte  flow). 
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7. 2 Plane  Poiscuille  Flow 

The  linear  stability  theory  was  applied  to  a plane  Poiseuille  flow 
of  centerline  velocity  Um  and  gap  width  25.  The  fluid  was  assumed  to 
obey  the  Oldroyd  constitutive  equation 


S1*  + 


1 


6S 

5t 


= 2p  (E1-^  + X 


SE 


2 5t 


-) 


(7.2) 


where  E1^  is  the  strain  rate  tensor,  A is  the  stress  relaxation  tine, 

1 

X2(<  A^)  is  the  strain  rate  relaxation  time,  and  o/5t  is  the  Oldroyd 

convective  derivative.  The  stability  equations  were  numerically 

integrated  by  modifying  the  existing  computer  technique  for  a Newtonian 

fluid.  The  resulting  neutral  stability  curves  (see  Fig.  4.8,  in  which 

A = A U /S  and  B = A U /6  are  two  elastic  parameters)  indicate  that 
1 m 2 m 

(i)  The  presence  of  viscoelasticity  destabilizes  the  flow. 

(ii)  At  a fixed  A , an  increase  of  A (<  A ) stabilizes  the  flow. 

1 2 1 

This  of  course  is  expected  since  the  Newtonian  results  should  be 
approached  as  A^  -*■  A^. 

(iii)  For  small  A and  A it  is  only  (A  -A  ) which  is  important, 

12  12 

not  A^  and  A separately. 

(iv)  The  critical  wave  number  increases  (that  is,  the  most 
dangerous  waves  become  shorter)  as  the  fluid  elasticity  is  increased. 

(v)  Our  results  do  not  agree  with  the  results  of  Tlapa  5 Bernstein 
(1970)  for  the  Maxwell  fluid  (A  = 0) ; an  examination  of  their  computer 
program  revealed  at  least  two  errors. 

It  has  been  noted  that  the  results  of  al 1 the  investigators,  who 
assumed  different  constitutive  equations  with  different  normal  stress 
behavior  and  Weissenberg  effect,  lead  to  an  identical  conclusion  for  the 
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plane  Poiseuille  flow.  This  fact  has  been  shown  to  be  a result  of  two- 
dimensionality  of  the  problem. 

7.3  Circular  Couette  Flow 

Here  again  the  linear  stability  theory  was  applied  to  the  flow  in 

the  annular  region  between  two  circular  cylinders,  with  and  R2  as  the 

radii  of  the  inner  and  outer  cylinders  respectively,  and  as  their 

angular  velocities,  and  5 = R -R^.  The  small  gap  approximation  was 

used,  and  the  fluid  was  assumed  to  obey  the  Oldroyd  constitutive 

equation  (7.2).  The  stability  equations  were  solved  by  means  of  an 

expansion  in  terms  of  a suitable  set  of  orthonormal  functions.  Numerical 

values  were  computed  for  two  cases,  namely  (i)  the  cylinders  rotating  in 

the  same  direction  with  the  same  angular  velocity,  that  is,  a = P.^/n^-1  = 

0,  and  (ii)  outer  cylinder  stationary,  that  is,  a = -1.  The  critical 

2 3 2 2 

values  of  the  nondimensional  Taylor  number,  T = -2a £1  5 R p /y  , were 

1 1 

found  by  drawing  the  neutral  stability  curves.  The  final  results  are 
given  as  curves  of  the  critical  Taylor  number  T£  vs.  A,  for  various 
values  of  B (Fig.  5. 4-5. 5),  and  the  critical  wavenumber  kc  vs.  A,  for 
various  values  of  B (Fig.  5. 6-5. 7). 

The  results  indicate  that  the  conclusions  (i)  through  (iv)  of  the 
last  section  are  also  valid  for  this  flow. 

The  analytical  results  cannot  be  compared  with  the  existing  experi- 


1 


mental  results  on  the  stability  of  circular  Couette  flow  of  polymer 
solutions,  since  the  experimental  results  of  different  observers  are 
sharply  contradictory. 
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7.4  Comparison  of  the  Three  Stability  Problems  Studied 

For  many  different  constitutive  equations  the  plane  Couette  flow 
has  predicted  a stabilization  due  to  elasticity,  while  the  plane 
Poiseuille  flow  has  predicted  a destabilization.  The  insensitivity  to 
the  constitutive  equation  has  been  shown  to  be  a result  of  two- 
dimensionality.  However,  it  is  not  easy  to  explain  why  the  plane  Couette 
flow  should  be  stabilized  and  the  plane  Poiseuille  flow  should  be 
destabilized.  The  obscurity  of  the  stability  behavior  of  these  two 
flows  are  well-known  for  the  Newtonian  fluid,  for  which  even  the 
mechanisms  responsible  for  causing  the  instability  are  little  understood. 
It  is  quite  possible  that  these  mechanisms,  whatever  they  are,  are 
modified  by  elasticity  differently  for  the  two  flows.  And  of  course 
there  is  the  experimental  evidence  of  the  Toms  phenomenon  that  the  plane 
Couette  flow  should  be  stabilized,  and  the  conclusion  of  destabilization 
of  the  plane  Poiseuille  flow  agrees  with  all  the  other  theoretical 
investigations  of  the  same  flow. 

In  presence  of  three-dimensional  disturbances  there  is  no  proof  of 
the  insensitivity  to  the  constitutive  equations,  and  the  results  of  the 
circular  Couette  flow  show  that  the  conclusions  are  indeed  dependent  on 
the  constitutive  equations  used.  All  the  investigations  (including  the 
present  one)  which  used  a constitutive  equation  predicting  a positive 
free  surface  Weissenberg  effect  predict  destabilization,  while  the  one 
investigation  of  Fong  (1965)  which  used  a constitutive  equation  having 
a negative  free  surface  Keissenberg  effect  predicts  stabilization. 

Since  a real  viscoelastic  fluid  generally  displays  the  positive  effect, 
the  conclusion  of  destabilization  seems  to  be  valid. 
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7.5  Normal.  Stresses  and  Weissenberg  Effect 

It  has  been  found  that  the  sign,  and  the  magnitude,  of  the 
Weissenberg  effect  may  depend  on  the  method  of  observation,  namely 
through  the  shape  of  the  free  surface  or  through  the  height  difference 
of  radial  tappings. 

The  Oldroyd  A fluid  and  the  Walters  A'  fluid  (both  long  and  short 
memories)  display  the  negative  Weissenberg  effect  as  far  as  the  free 
surface  is  concerned,  but  they  would  display  the  opposite  effect  when 
observed  through  radial  tappings.  The  Oldroyd  B as  well  as  the  Walters 

B1  fluid  (both  long  and  short  memories)  display  identical  positive 

* 

effect  in  both  types  of  observations.  The  second-order  fluid  would 
display  the  positive  Weissenberg  effect  when  observed  through  radial 
tappings,  but  it  may  display  either  effect  when  the  free  surface  is 
observed.  A simple  experiment,  namely  that  of  observing  the  Weissenberg 
effect  by  the  two  methods,  is  suggested  for  measuring  the  normal  stress 
functions  of  a viscoelastic  liquid. 

7.6  Epergy  Consideration  of  Constitutive  Equations 

It  has  been  mentioned  in  Section  7.1  that  in  order  for  the  energy 
method  of  Orr  to  become  applicable  to  a second-order  fluid,  it  was 
necessary  to  bring  in  the  concept  of  the  elastic  potential  energy,  and 
to  discover  an  expression  for  it. 

For  a purely  viscous  liquid  all  the  work  spent  in  deforming  a 
material  particle  is  dissipated  into  internal  thermal  energy.  It  is 
expected  that  in  a viscoelastic  liquid  a part  of  the  shearing  deforma- 
tion work  must  be  stored  as  elastic  energy.  Consequently  the  rate  of 
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deformation  work  in  a viscoelastic  liquid  has  been  written  as  a sum  of 
dissipation  and  a rate  of  increase  of  elastic  potential  energy: 

ViJ  ■ T * TT1  <7-3> 

where  T\  ^ is  the  stress  tensor,  Lh  is  the  velocity,  Y is  the  dissipation, 
and  PE  is  the  elastic  potential  energy.  Expressions  have  been  found  for 
Y and  PE  of  the  various  commonly  used  constitutive  equations,  namely  the 
second-order  fluid,  Oldroyd  fluid,  Maxwell  fluid,  and  Walters  fluid  with 
short  memory. 

It  has  been  found  that  the  dissipation  Y for  the  Oldroyd  and  Maxwell 
fluids  is  always  positive,  but  the  dissipation  for  the  second-order  and 
the  Walters  fluids  is  not  necessarily  so.  However,  it  has  also  been 
found  that  Y is  positive  for  the  latter  two  fluids  for  all  two- 
dimensional  motions,  and  also  for  all  slow  flows;  this  conclusion  is 
consistent  with  the  fact  that  these  two  constitutive  equations  were 
derived  as  slow  flow  approximations  to  more  general  constitutive 
equations . 

7.7  Some  Unsteady  Parallel  Flows 

The  following  three  parallel  flow  problems  of  the  Oldroyd  fluid 
have  been  solved: 

(i)  Flow  due  to  a single  flat  plate  oscillating  in  an  infinite 
fluid. 

(ii)  Flow  in  the  space  between  two  parallel  plates,  one  of  which 


is  oscillating. 
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i(kx2-o)t) 

(iii)  A shear  wave  motion  given  by  = Ae  , with  k real 


and  w = ur  + iw^. 

No  startling  facts  have  been  found  in  the  first  two  problems. 
However,  an  interesting  fact  was  found  for  the  third  problem  for  which 
it  was  found  that 


+ f\, 


oi  = 


i(l+vk  X2)  + /4AjV  2 - (l+k2vA2)2 

__ 


(7.4) 


where  v is  the  kinematic  viscosity.  Two  cases  are  therefore  possible, 

2 2 2 
according  as  4A^k  is  more  or  less  than  (1  + A2vk  ) . In  the  first 

case  (large  elasticity)  > 0 and  ok  < 0,  so  that  the  waves  are  propa- 
gated as  well  as  damped  in  time.  In  the  second  case  (small  elasticity) 
o>r  = 0 and  oi^  < 0,  so  that  the  waves  do  not  propagate  but  are  "standing" 
and  damped,  as  in  Newtonian  fluids. 

For  a Maxwell  fluid,  the  above  phenomenon  is  explained  as  follows: 
in  one  case,  the  propagation  speed  of  elastic  waves  is  greater  than  that 
of  the  viscous  diffusion,  so  that  there  is  a net  propagation;  in  the 
other  case,  the  viscous  diffusion  overcomes  the  elastic  wave  speed,  so 
that  there  is  no  net  propagation. 


7.8  Final  Remarks 

The  soundness  of  the  basis  of  the  entire  field  of  viscoelastic  flows 
seems  questionable  to  us.  For  example,  the  conclusions  sometimes  depend 
on  the  constitutive  equations  used.  Leaving  aside  subtle  questions  like 
the  hydrodynamic  stability,  they  display  different  types  of  normal 
stresses  and  IVcissenberg  effects  in  the  simplest  of  all  flows — a steady 
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viscometric  flow.  It  is  highly  desirable  that  the  normal  stress 
functions  of  a real  viscoelastic  liquid  be  thoroughly  investigated 
experimentally.  The  various  existing  studies  on  this  subject  have 
failed  to  agree  even  on  the  sign  of  these  functions.  It  is  suggested 
that  many  more  studies  should  be  undertaken  to  settle  this  question. 

It  may  be  possible  that  different  materials  have  different 
behavior.  Some  constitutive  equations  may  then  be  appropriate  for  one 
fluid  and  others  for  another.  From  a thorough  experimental  study  of 
the  normal  stress  functions  of  different  materials,  a rational  ground 
regarding  the  choice  of  the  appropriate  constitutive  equation  for  that 
particular  material  would  emerge. 

Lastly,  too  much  confidence  should  not  be  placed  in  our  explanation 
of  the  thickening  of  the  viscous  sublayer  during  the  Toms  phenomenon,  as 
given  in  Sections  3.7  and  7.1.  Firstly,  the  disturbances  in  the  sub- 
layer are  observed  to  be  of  the  type  of  streamwise  vortices,  for  which 
we  have  demonstrated  a destabilization  of  the  plane  Couette  flow  (and 
therefore  a thinning  of  the  viscous  sublayer) , if  the  (second-order) 
fluid  i?  to  display  the  positive  Weisservberg  effect.  Secondly,  there 
are  some  deficiencies  in  Taylor's  argument,  as  pointed  out  by  Lumley  5 
McMahon  (1963).  The  author  therefore  feels  that  there  may  not  be  any 
connection  between  the  Toms  phenomenon  and  the  stability  of  visco- 
elastic flows;  the  explanation  of  the  phenomenon  may  lie  elsewhere. 
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Three  problems  in  the  stability  of  viscoelastic  flows  have  been  theoreti- < 
cally  investigated. 

The  case  of  plane  Couette  flow  has  been  solved  by  the  classical  energy 
method  of  Orr.  \For  the  method  to  be  applicable  it  has  been  found  that  the 
concept  of  elastic\ootential  energy  of  a viscoelastic  liquid  has  to  be  intro-  j 
duced,  and  an  expression" ha's  "been  found  for  it.  * For  two-dimensional  dis- 
turbances of  any  magnitude  the  presence  of  elasticity  has  been  found  to 
stabilize  the  flow.  This  result  suggests  that  the  viscous  sublayer  thickness 
must  increase  during  the  so-called  Toms  phenomenon.  It  has  been  found  that  the 
results  are  identical  for  three  fluids,  namely  the  second-order  fluid  and  the 
Walters  fluids  A'  and  3',  although  the  normal  stress  behavior  of  these  fluids 
are  quite  different.  This  insensitivity  to  the  constitutive  equation  has  been 
shown  to  be  a result  of  the  two-dimensionality  of  the  problem. 

Linear  stability  theory  has  been  used  to  investigate  the  stability  of  two 
flows  of  the  Oldroyd  fluid,  namely  the  plane  Poiseuille  flow  and  the  circular 
Couette  flow.  For  both  flows  the  presence  of  elasticity  has  been  found  to  de- 
stabilize the  flow  and  to  increase  the  critical  wavenumber 

As  an  aside,  it  has  also  been  shown  that  the  sign,  and  magnitude,  of  the 
Weissenberg  effect  may  depend  on  the  method  of  observation,  namely  through  the 
shape  of  the  free  surface  or  through  the  height  difference  of  radial  tappings.  \ 
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